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Section 1 
INTRODUCTION 

This report summarizes work performed during the second phase of an effort to 
develop a computer- implemented numerical method for predicting the flow charac- 
teristics and performance of three-dimensional jet engine exhaust nozzles. The 
objective of developing a method for computing the internal and external viscous 
flowfield of an isolated nozzle has been met. The approach is based on using 
an implicit numerical method to solve the unsteady Navier-Stokes equations in 
a boundary-conforming curvilinear coordinate system to obtain the desired time- 
asymptotic steady state solution. Flow turbulence effects are simulated by 
means of algebraic turbulence models for the effective turbulent eddy viscosity 
and Prandtl number. A detailed description of the equations and boundary con- 
ditions and of the numerical method have been presented in an earlier report 
[Ref. 1], along with a general discussion of turbulence models appropriate to 
the various sub-regions of the nozzle flowfield. 

The present final report describes work performed since Reference 1 was written. 
Recent modifications and improvements to the original numerical algorithm are 
presented in Section 2. Section 3 gives the equations that are used to derive 
the nozzle performance parameters such as thrust and discharge coefficient 
from the computed flowfield data. The final formulation of the turbulence 
models that are used to simulate flow turbulence effects is presented in 
Section 4. Section 5 presents the results of numerical experiments performed 
to explore the effect that various parameters in the numerical method have on 
both the rate of convergence to steady state and on the final flowfield solu- 
tion. Detailed flowfield predictions for several three-dimensional nozzle con- 
figurations are presented in Section 6 and compared with experimental wind 
tunnel data. 

The numerical method is embodied in a set of three computer codes: RGRIDD, 

NOZLIC, and N0ZL3D. The RGRIDD code constructs the curvilinear coordinate 



system and computational grid numerically for nozzles of complex geometric 
configuration. The NOZLIC code generates a set of flowfield initial condi- 
tions on this grid that are used to start a flow computation. The N0ZL3D 
code performs the actual flowfield computation and evaluates the nozzle per- 
formance characteristics. A user's guide to the operation of these three 
codes is contained in a separate volume [Ref. 8]. 

In the sections that follow, all equations are cast in dimensionless form. 
Distances are referred to a reference length, velocities are referred to the 
speed of sound at some reference state, viscosity y is referred to the 
molecular viscosity at the reference state, and individual state variables 
such as density p , pressure p , and temperature T are referred to their 
values at the reference state. The reference state is that at which the 
Reynolds number Re is defined in terms of the reference length using the 
reference sound speed as the characteristic velocity. 
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Section 2 

MODIFICATIONS TO THE NUMERICAL METHOD 
2.1 BRIEF OUTLINE OF THE ORIGINAL FORMULATION 

The original formulation of the numerical method [Ref. 1] starts with the 
strong conservation-law form of unsteady Navier-Stokes equations in a Cartesian 
base coordinate system (xyz). The equations are transformed to a boundary- 
conforming curvilinear coordinate system U,n>c), and take the non-dimensional 
form 


q T + f 5 + 9 n + - Re' 1 (e n + ) ( 2 . 1 ) 
where q = Jq" (2.2) 
J is the Jacobian of the inverse transformation 


J = 9(x,y,z)/a(£, nj c) 


(2.3) 


q is a vector of conserved variables whose components are the density, the 
three Cartesian components of the momentum flux vector p V, and the total 
energy per unit volume, f, g, and h are inviscid flux vectors; and the terms 
on the R.H.S. of Eq. (2.1) that are inversely proportional to the Reynolds 
number Re represent viscous transport processes. Each of the inviscid flux 
vectors is a linear combination of the flux vectors ?, g, ft associated with 
the Cartesian coordinate direction x, y, z, respectively. The coefficients 
of these linear combinations are the metrics of the coordinate transformations 
e(x,y,z,t), n(x,y,z,t), e(x,y,z,t). For example. 


f = E t q + « x f + ? v g + C, & 


E - J £ = yz-yz 

x s x •'n ? n 


h = J «y = z n " z c x n 


etc. 
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The non-dimensional flow variable vector q and the inviscid flux vectors 

■ 4 - 

f, g, h are defined as follows in terms of the Cartesian components u, v, w, 
the density p, the pressure p, and the total energy per unit volume e: 


q = (p, pu, pv, pw, e) T 

f = [pu, p'+ pu , puv , pUW, u( p 1 + e)] 

g = [pv, pvu, p'+ pv 2 , pvw, v( p 1 + e)] T 

h = [pw, pwu, pwv, p'+ pw , w(p' + e) ] 


where p' = p/y 


and y is the specific heat ratio. 


Similarly, each of the viscous terms 6 , w is a linear combination of the 
viscous terms associated with each of the Cartesian coordinate directions. The 
reader is referred to Section 2 of Reference 1 for the mathematical equations 
that define the viscous terms. We note only that Eq. (2.1) represents the 
parabolic approximation to the Navier-Stokes equations wherein viscous terms 
associated with the streamwise coordinate E, are neglected. 


To solve Equation (2.1), the transformed space (£, n, t) is covered by a 
uniform grid (Cj, n ^, 5 ) such that peripheral grid points lie on the bound- 
aries of the space. The spatially differenced form of the equations is derived 
to second order accuracy in the mesh spacings A£ , An, A? by means of the 
Finite-Volume Method. To each interior grid point of the transformed space 
there corresponds a cell of volume At/ = A£ An A? that encloses the point. 

The difference equations that apply at the point are derived by integrating 
Eq. (2.1) over the cell volume. This leads to difference equations of the form 

q d n dc + Fj <5j f + <5 k g + h = ... (2.4a) 



At/ = A? An A? 


(2.4b) 
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where the viscous terms have been suppressed for brevity, and where u. 6. 

J J 

is the centered spatial differential operator for the Sj direction 




(2.5) 


and the central difference operators for the k and a directions are defined 
similarly. 

Within the second order spatial accuracy of the remaining terms, the volume 
integral that appears in the first term of Eq. (2.4) may be represented in 
terms of the cell -averaged value 


^ ( A ^ -1 f ^ d? dri dc 


( 2 . 6 ) 


which is centered at the grid point jkJt itself since the grid point is located 
at the centroid of the cell that surrounds the point. 

The final implicit space-time difference equations that govern the change 

A <N n n*|-l tt 

Aq = q - q over a time step At = x - x are obtained by evaluating 

n+1 

the spatially differentiated terms at the advanced time t , performing a 
first-order Taylor series expansion about the solution at time x n , and 
factoring the implicit operator. This yields the implicit ADI sequence 

A u A A A ^ M 

(I + At jij 6^ F) Aq = -(yj ^ + y |< 6 k 9 + ^ At + **‘( 2 * 7a ) 

(I + At <5^ G + ...) n Aq = Aq (2.7b) 

(I + At <5^ H +...)^ Aq = Aq (2.7c) 

where F, G, H are Jacobian matrices 


f - , e-|a , h-£ 

3 q 3 q d q 
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and where the viscous terms again have been suppressed for brevity. Each step 
of the ADI sequence in Eq. (2.7) involves solving a block-tridiagonal linear 
system of equations to obtain the solution at interior grid points. 

For grid points located on the boundaries of the computational space C,n>C s 

one or more of the five scalar components of Eq. (2.1) are replaced by algebraic 

boundary conditions. If the latter are nonlinear in time, they are linearized 

_*n 

by a first-order Taylor series expansion about the solution q . This yields a 
linear algebraic subsystem of the form 

M n q = in 11 (2,8) 

The remaining scalar components of Eq. (2.1) (i.e,, those scalar flow equations 
that have not been replaced by the aforementioned algebraic boundary conditions) 
are differenced by applying the finite- volume method in the same way as des- 
cribed above for interior points. However, to a typical boundary point such 
as j = 1 , there corresponds not a full cell, but rather a half-cell whose 
width in the coordinate direction normal to the boundary is only half the width 
AC of an interior cell. The counterpart of Equations (2.4) take the form 

37 f q dc dn dc + Aj f + 6 ^ g + y £ h = . . . s j = 1 (2.9a) 

IV 


AV = (ac/ 2) An AC 


(2.9b) 


where Aj is the forward difference operator such that Aj = fj + -j - fj. The 
volume integral in the first term of Eq. (2.9a) still may be represented in 
terms of the value q* at the cell centroid 


q* = (At/) 


'l 


q dc dn dc 


( 2 . 10 ) 


but the centroid is not located at the boundary grid point j = 1. However, 
within the spatial order of accuracy of Eq. (2.9), the value q* may be 
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evaluated by linear interpolation between the boundary point j = 1 and the 
adjacent interior point j = 2. 

q* = (I + i Aj) q j = 1 (2.11) 

where I is the identity operator. 

Equation (2.9a) is time-differenced implicitly and linearized in the same way 
as is Eq. (2.4) for interior points. The appropriate scalar components of the 
resulting equation are replaced by the linearized algebraic boundary conditions 
(2.8) and the implicit operator is factored to obtain the counterpart of the 
sequence (2.7) that applies at boundary points. 


2.2 MODIFIED FORMULATION OF THE IMPLICIT ALGORITHM 

The original formulation of the algorithm as presented in Reference 1 and 
summarized in the preceding subsection has two deficiencies: (1) the algorithm 

is valid only if the curvilinear coordinate transformation has no singular 
points where the Jacobian J vanishes, and (2) unacceptably large truncation 
errors can arise at grid points situated along the lines of intersection be- 
tween boundary surfaces of the computational domain. 

As an example where a singular transformation arises, consider the internal 
flow in an axi symmetri c nozzle whose axis of symmetry coincides with the 
Cartesian x axis, and whose interior wall is of radius r = r w (x), 

0 5 x £ L. The first quadrant of the flow region interior to the nozzle can 
be mapped onto a rectangular parallelepiped in a right-handed curvilinear 
coordinate system ti, ? by the transformation 


5 = 

€(x) , 

0 £ x £ L 


n = 

-e 

0 < 0 < ir/2 

(2.12) 

S = 

s(r) , 

0 5 r 1 r w (x) 
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where 


e = tan”^ y/z 
r = 

One can verify easily that the Jacobian J of this transformation vanishes 
at the axis of symmetry y = z = 0. The axis of symmetry maps onto the face 
S - 0 of the parallelepiped in the transformed space £, n» hence, the 
transformation is singular at each point of that face. 

The numerical mapping technique presented in Section 3 of Reference 1 also can 
generate curvilinear coordinate systems that have isolated singularities. For 
example, the quasi-elliptical mapping depicted in Figs. 3-3 and 3-4 of 
Reference 1 has a singularity at the point A that corresponds to the focus 
of the ellipse. The quasi-rectangular mapping given in Fig. 3-6 of Ref. 1 
for the flow regions interior and exterior to a nozzle whose cross-section is 
super-elliptical also has singularities at the points of maximum curvature of 
the internal and external super-ellipses. 

The numerical algorithm can be modified to handle isolated singularities such 
as those in the foregoing examples. Two modifications are necessary: first, 

as dictated by the finite-volume method, the Jacobian J at any grid point that 
coincides with a coordinate singularity must be computed as a cell-averaged 
quantity; and, second, the ADI sequence (2.7) must be modified to yield Aq 
directly, rather than Aq, which involves the Jacobian as a factor. We shall 
deal with these modifications in order. 

In general, the flow variable vector q is regular and non-zero even in the 
neighborhood of a coordinate singularity, whereas the compound quantity q in 
in Eq. (2.2) vanishes along with J at the singularity itself. Thus, Eq. (2.6) 
is a valid representation of the volume integral that appears in the first term 
of Eq. (2.4a) only at grid points where the 5, n» t coordinate transformation 
is non-singular. At a singular point, the quantity q vanishes locally, 
whereas the volume integral is non- zero because it includes contributions from 
all regular points within the finite-volume cell that encloses the singular 
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point. Since q itself is regular, the volume integral in Eq. (2.4a) can be 
evaluated by applying the mean value theorem 


( ai /)- 1 



q d£ dn d? 



(2.13) 


where the cell-averaged Jacobian J* is given by 


J* * (AlO 


•'A l / 


J d£ dn dc 


(2.14) 


The latter is always non-zero for any cell of finite volume At/. If the singu- 
larity coincides with an interior grid point, then the integral in Eq. (2.14) 
can be evaluated analytically to second-order accuracy by introducing a local 
Taylor series expansion for the function J( e, n» c). If the singularity 
coincides with a boundary grid point, say j = 1, then the counterpart of 
Eq. (2.13) is 

(AtO' 1 f q d? dn d; - q* J* (2.15) 

{v 

where both q* and J* represent cel 1 -averaged values at the cell centroid, 
and can be computed individually by linear interpolation between the boundary 
point and the adjacent interior point (cf. Eq. (2.11)). 

The alterations necessary to permit the direct computation of Aq from the 
algorithm are merely a special case of the alterations that permit the use of 
grids that move as a function of time (Ref. 2). One need only expand the term 

'S 

Aq as follows 


Aq 


->n 

q 


AJ + J 


n+1 


Aq 


and redefine the Jacobian matrices in Eq. (2.7d) as 
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F = 


9f 


= aa 


ah 

3 ? 


Upon factoring the implicit operator in the implicitly time-differenced and 
linearized version of Eq. (2.4a), one obtains the ADI sequence 


(J n+1 I + At Uj F n )Aq = -q n Aj -(y^ f + y k 6 ^ g 


(2.16a) 


(J 

(J 


n+1 

n+1 



+U^ 5 

z h) Ax + ... 


I + At y^ 6^ + 

...) a3 * 

- I n + 1 A^ 

= J Aq 

(2.16b) 

I + At y^ <5^ H* 1 + 

...) a5 

,n+l 

= J Aq 

(2.16c) 


where AJ = 0 and J n+1 = J n if the grid is stationary in time, i.e., if the 
transformation (x,y,z)+(£;, q, t;) is independent of time. 


It is important to note that the form given in Eq. (2.8) must be retained when 
applying algebraic boundary conditions at boundary grid points, or the values 
of Aq and of Aq at boundary points will become inconsistent with those 
at interior points where Eqs. (2.16) are employed. That is, for a stationary 
grid, the boundary conditions must be written as 

J M n Aq = in n (2.17) 


The modified algorithm given in Eqs. (2.16) and (2.17) should be valid at grid 
points that coincide with coordinate singularities as well as at regular points. 
However, numerical experiments for axi symmetric flow yield poor numerical re- 
sults for the flow variables at grid points situated along the singular axis 
of symmetry, although the solution is accurate at all other points. We con- 
jecture that the poor results at the symmetry axis result from a locally large 
truncation error, inasmuch as care was taken to incorporate the symmetry pro- 
perties of the flow variables and of the Cartesian coordinates into the computa- 
tion of the metrics and of the explicit fourth order smoothing terms in the 
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neighborhood of the axis. To obtain an accurate solution at points on the 
axis, we have found it necessary to employ a time-lagging approach in which 
the flow variables at axis points are extrapolated from those at adjacent 
points following each time step. Even functions such as the temperature, 
density, and axial velocity component are extrapolated from those at the two 
adjacent interior points using a second degree polynomial whose gradient is 
zero at the axis, and odd functions such as the transverse velocity components 
are set to zero at the axis. 

Locally large truncation errors also are incurred at other types of coordinate 
singularity, as well as exceptional but non-singular points of the curvilinear 
coordinate system. Examples of the latter are grid points that lie along the 
curves of intersection between two families of coordinate surface n = const, 
and c = const, that represent nozzle wall boundaries, such as the interior or 
exterior corner of a so-called "two-dimensional 11 nozzle that has flat walls 
and a rectangular cross-section. Such a corner point has the characteristic 
feature that the finite-volume cell associated with the point is either a 
quarter-cell (interior corner) or a three-quarter-cell (exterior corner). To 
avoid large numerical inaccuracies at such exceptional points, we have employed 
a time-lagging approach similar to that outlined above for singular points of 
an axi symmetric flow. For example, at axial corners where two walls intersect, 
we merely interpolate the temperature and density from the nearest neighboring 
wall points at the end of each time step. 


* When the Cartesian x-axis coincides with the flowfield axis of symmetry, 
the coordinates y, z and their associated velocity components, v, w 
are odd functions of position relative to the axis, whereas x and all 
other flow variables are even functions. 
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2.3 MODIFIED SCHEME FOR SUBSONIC INFLOW AND OUTFLOW BOUNDARY CONDITIONS 


2.3.1 Inflow Boundary Points 

The scheme outlined in Section 4.2.2 of Reference 1 has been modified slightly 
to avoid an obvious inconsistency at inflow boundary grid points that lie on 
the nozzle walls. In the original inflow boundary scheme, four algebraic 
boundary conditions are applied that specify the total pressure, the total 
temperature, and the two direction cosines of the velocity vector. The fifth 
relation that is required to close the system of equations governing the five 
unknown components of the flow variable vector q is obtained from an implicit 
finite-volume discretization of the mass conservation equation (the first 
scalar component of Eq. (2.1)). This use of the mass conservation equation 
to determine the density is invalid at wall points because it is inconsistent 
with the density that is implied by the total pressure and total temperature 
through the equation of state. This follows from the fact that total and 
static temperatures are equal and total and static pressure are equal at wall 
points where the velocity vanishes. Thus, the algebraic boundary conditions 
alone are sufficient to determine the flow variables at wall points of the 
inflow boundary. Note, however, that when the temperature T w specified as 
a boundary condition at nozzle walls, the inflow stagnation temperature boundary 
condition must be equal to T w at points where the walls intersect the inflow 
boundary. 

2.3.2 Outflow Boundary Points 

For cases where one is interested in computing only the flow internal to the 
nozzle, the outflow boundary is positioned at the nozzle exit plane. When the 
ambient pressure outside the nozzle is sufficiently high relative to the internal 
flow stagnation chamber pressure, the flow will be wholly subsonic at the exit 
plane, and the ambient pressure must be imposed as a boundary condition on the 
static pressure at the outflow boundary {Ref. 1, Section 2.4.6]. In the original 
implicit algorithm for outflow boundary grid points, this boundary condition is 
linearized and used in place of the u-momentum equation. However, this results 
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in an inconsistent set of equations for the flow variables at grid points 
located at the nozzle wall, where additional algebraic boundary conditions 
on the velocity components and on either the wall temperature or heat flux 
also are imposed. The inconsistency becomes apparent when one observes that, 
for an adiabatic wall, the continuity and energy equations determine the gas 
density and temperature at the wall. The corresponding wall pressure implied 
by the equation of state in general will be inconsistent with the boundary 
condition on the exit plane pressure. A similar inconsistency exists at wall 
points when the wall temperature is specified as a boundary condition. The 
inconsistency can be avoided by retaining the u-momentum equation at all out- 
flow boundary points, and using the imposed exit pressure boundary condition 
in place of the energy conservation equation. In addition, when the wall 
temperature T w is specified as a boundary condition, the algebraic equation 
T = T w is used in place of the continuity equation at wall points of the out- 
flow boundary. The physical justification for this procedure is that the gas 
density at the wall is determined completely by the exit pressure and wall 
temperature boundary conditions alone. 


2.4 IMPLICIT DISSIPATION 

The implicit algorithm permits the use of large time steps At without numeri- 
cal instability, and makes it possible to attain the steady-state solution of 
Eq. (2.1) in fewer time steps than would an explicit algorithm. However, 
hundreds of time steps usually are required to achieve convergence to steady 
state. In an effort to speed convergence, we have introduced artificial 
dissipative terms into the implicit one-dimensional operators on the L.H.S. 
of 'Eqs. (2.16). These dissipative terms are similar to those employed by 
Steger [Ref. 3], except that they are differenced in a conservative fashion 
and obey homogeneous boundary conditions. This ensures that the dissipative 
terms do not alter the global conservation properties of the difference 
equations [Ref. 1, Section 4.4], 

The form of the dissipative term for the j coordinate direction is 
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d 


, j = 1 

, 1 < j < j, 


( 2 . 18 ) 


a Ej J 6 . (Aq) 

= {a J (Aq) 
i-a ET 1 J 6 . (Aq) 


max 


J = J 


max 


where a is a constant; E is the classical shift operator which is defined 

+m 

such that for any mesh function f., f = f . ; 6 is the classical central 

i Ji J J±m 

difference operator 6.f = (E? - ET 2 ) f; and J is the transformation 

vJ J J 

Jacobian in Eq. ( 2 . 3 ). The central member of Eq. (2.18) applies at interior 
points, and the first and last members apply at boundary points that are not 
situated at flowfield symmetry planes. For such symmetry boundary points, the 
dissipative terms are modified to account for the symmetry properties of q 
and of J. With the addition of the dissipative term (2.18), the implicit 
operator on the left side of Eq. (2.16a) for interior points assumes the form 

(J n+ ^ I + a 5j J n+ ^ Sj I + At y. <5 j F n ) Aq = ... 

Similar dissipative terms are added to the implicit operators in Eqs. (2.16b) 
and (2.16c) for the n and c coordinate directions. These dissipative terms 
do not affect the final steady state solution, because Aq vanishes at steady 
state. Furthermore, they do not alter the unconditional numerical stability 
of the algorithm. 


The introduction of the implicit dissipative terms is equivalent to appending 
terms of the form 


a(A£)^ 



(2.19) 


to the right side of the Navier-Stokes Equations (2.1). That is, a Taylor 
series expansion of Eq. (2.16) with the dissipative terms yields a modified 
partial differential equation which, to lowest order in At, a?. An, Ac is 
identical to Eq. (2.1) with three additional terms of the form (2.19), one 
for each coordinate direction. The steady state solution of the Navier-Stokes 
equations is unaffected by the additional terms because all time-derivative 
terms vanish at steady state. 
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Numerical experiments indicate that the artificial implicit dissipative terms 
have a favorable effect on the convergence rate when dissipation coefficient 
a is of the order of unity (see Section 5.4). However, initial experimenta- 
tion with the dissipative terms displayed a tendency to produce oscillations in 
the spatial distribution of computed flow variables across a subsonic inflow 
boundary. To damp these oscillations, which have a wavelength of twice the 
mesh spacings, An, A?, we have found it necessary to filter the computed 
flow variable distribution at the inflow plane following the second and third 
steps of the ADI sequence (2.16) This is accomplished by the one-dimensional 
low-pass filters 


q -h 1 

Aq = y 2 Aq 
Aq = y^Aq 


( 2 . 20 ) 


where y denotes the classical central averaging operator 



Eqs. (2.20) apply only at interior points of the inflow plane. No filter is 
applied at the boundaries of the inflow plane unless those boundaries coincide 
with flowfield symmetry planes. For example, the intermediate solution Aq is 
not filtered over k at boundary points such as k = 1, k max unless those points 
are located at symmetry planes, in which case the appropriate filter is 
Aq = 2 P k Aq , where the plus sign applies at k = 1 and the minus sign applies 
at k ^max' 


The described filters are applied only to the first, second, and fifth components 
of the flow variable vectors Aq and Aq . The transverse momentum components 
Apv, Apw are recomputed from the filtered component Apu using the inflow 
boundary conditions on the direction cosines of the velocity vector (see Section 
2.3.1 above). 
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2.5 IMPROVED NUMERICAL GRID GENERATION TECHNIQUE 


A general technique for generating a boundary-conforming curvilinear coordinate 
system £,n,C and computational grid suitable for geometrically complex nozzle 
configurations has been given in Section 3 of Reference 1. In this technique, 
a three-dimensional grid is built up by constructing a sequence of two-dimensional 
grids in successive cross-sectional planes £(x) = const. Within each such y-z 
cross-sectional plane, the transverse coordinate system and computational 
grid is generated numerically as the solution to an elliptic boundary value 
problem governed by the following elliptic system of equations 


where 


“ ( Yrn +t Wn ) - 2 By r , c +-r' ( y n r +i'y ? ) = 0 

(2.21a) 

a ^n + % ) ' 2ez nC + ' r(z «' h| ' z 5 ) = 0 

(2.21b) 

a = y^+zj; 

(2.22a) 

6 = 

(2.22b) 

y = y 2+z ^ 

1 J T\ n 

(2.22c) 


Equations (2.21) are solved numerically on a uniform, rectangular grid n k > 
\ = (k-l)An ,k = 1,2,. . . 

= U-1)A? ,£, - 1,2,. . . 


to obtain the Cartesian coordinates (y,z) of the grid point in physical space 
that corresponds to each point in the computational space. The boundary 

values for Eq's. (2.21) are the y,z coordinates of grid points on the boundaries 
of the flow region in the physical domain. These boundary grid points may be 
distributed unequally along the boundaries in any fashion. The parameters 
ip,$ in Eq's. (2.21) are evaluated locally at the boundaries in terms of the 
given boundary values (y,z) by using limiting forms of the elliptic equations. 
These parameters then are interpolated into the interior of the domain from the 
boundaries, and the elliptic system is solved numerically by a standard 
successive line over-relaxation technique. This results in a grid point 
distribution throughout the physical domain that is controlled entirely by the 
priori selection of the grid point distribution along the boundaries of that 
domain. 
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The described procedure always yields a boundary-conforming transformation 
in which the boundaries of the flow region in the physical domain are coordinate 
curves ri = const, or c = const, of the curvilinear coordinate system. In the 
original technique (Ref. 1, Section 3), the parameters and were evaluated 
from the boundary values using limiting forms of Eq's. (2.21) that were obtained 
by assuming that partial derivatives with respect to the curvilinear coordinate 
transverse to the boundary vanish locally at the boundary. For example, for a 
boundary £= const., the ^-differentiated terms were dropped from Eq's. (2.21) 
to obtain the limiting forms 

= 0 (2.23a) 

z nn + ‘ pz r] = 0 < 2 - 23t> ) 

These equations then were used to evaluate the parameter cp locally at each grid 
point on the boundary in terms of the boundary values of y,z by replacing the 
differential operators by central difference operators 


v ~ u k 6 k y (2.24a) 

An 

y ^ 6 k y (2.24b) 

To avoid locally large numerical errors at points where [y^ is small, Eq. (2.23a) 
was used only at points where Jy r] J^jz r) J, whereas Eq. (2.23b) was used at points 
where | y^ j c | z n | (Ref. 1, Section 3). 

The described procedure for evaluating the parameters from the boundary values 
was found to yield excellent computational grids for a variety of nozzle 
conf igurations. However, that procedure rests on a weak assumption, namely, that 
the derivatives in the direction transverse to the boundary can be dropped from 
the equations. One can show that these transverse derivatives actually vanish 
identically only when the boundary is both straight and is parallel to one of 
the Cartesian coordinate axes y or z. When this is not the case, the transverse 
derivatives are non-zero at the boundary. We recently have discovered that it is 
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possible to derive a universally valid limiting form of the elliptic system 
at the boundaries simply by imposing a local constraint on the angle of inter- 
section between the two families of coordinate curves £ = const, and n = 
constant. In particular one may impose the constraint that the two families 
be orthogonal everywhere along the boundaries*. This remarkable result can 
be proved as follows. 

Consider the case where we wish to evaluate the function v in Eq's. (2.21a,b) 
at a boundary £=£[ 3 = constant. Upon eliminating the function ip between the 
two equations, one obtains a single equation that can be cast in the form 

a [ z c^un + ^n)-yc( 2 nn +<pz n)] = " Y (^r) j ‘ 


Now, the ratio y^/z^ is merely the slope dy/dz of the family of coordinate curves 
n = const, that are transverse to the boundary curve c = C 5 . We are at liberty 
to impose the constraint that the transverse coordinate curves n = const, be 
locally straight (i.e. 3 have zero curvature) in the neighborhood of the boundary. 
This constraint may be stated in the form 



We now impose the further constraint that these transverse coordinate curves 
n = const, be locally orthogonal to the boundary c = £ 5 . The orthogonality 
condition may be found as follows. Let r = (y,z) denote the radius vector in 
the Cartesian y-z plane. Then the local tangent vector to a coordinate curve 
n = const, is 


* Note that this does not necessarily imply that the curvilinear coordinates 
will be orthogonal in the interior. 
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(2.27a) 


h = (y^) 

Similarly the local tangent vector to a coordinate curve £ = const, is 

r n = (y^) (2.27b) 

The two families of coordinate curves then are orthogonal if and only if 

r n -r 5 = 0 (2.28) 

The orthogonality condition (2.28) may be expressed in the form 

y n y c +z n z c =0 at ? = ? b (2.29) 

When we evaluate Eq. (2.25) at the boundary £ = the second term in brackets 
on the R.H.S. vanishes by virtue of Eq. (2.26); the first term in those 
brackets also vanishes since 3 = 0 by virtue of the orthogonality relation (2.29). 
The latter relation may be used to eliminate all ^-differentiated terms from 
the L.H.S. of Eq. (2.25). This yields a limiting form of the elliptic system 
that is valid at the boundary c = c b , an d that can be solved directly for the 
parameter <P 


V “ (y r |Y r in +z n z riTi ) on 4 = 4b (2.30) 

This represents a universally valid equation that can be used to compute the 
numerical value of <p at each grid point on the boundary in terms of the boundary 
values y,z once the differential operators are replaced by the difference 
oprators (2.24). The corresponding expression that determines the parameter 
along boundaries n = const, can be obtained directly from Eq. (2.30) by the 
substitution The values of the parameters throughout the interior of 

the n,c domain then are found by linear interpolation as in the original method. 
For example, <P(c,n) is computed from its values at the two boundaries c = const, 
by linear interpolation along lines n = const, in the rectangular computational 
domain n,£. This ensures that the final grid obtained from a numerical solution 
of the elliptic system (2.21) will reflect the boundary value distribution, and 
will have the desirable property that the two families of grid lines are locally 
orthogonal at the boundaries of the physical flow region. 
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It is instructive to explore the geometric interpretation of Eq. (2.30), which 
is used to evaluate the parameter <p along a boundary £ = Cb = const, in terms 
of the pre-assigned boundary values of y,z ; That equation can be re-cast in 
the form 


s nn +<ps n = 0 on c = Cb (2.31a) 

where s denotes the arc length along the boundary curve C = ?b 

ds = ^ dx^+dy 2 (2.31b) 

Eq. (2.31a) clearly possesses exponential solutions if <p is constant. Thus, 
the use of Eq. (2.30) to evaluate the parameter <P at each point along the boundary 
£ = ^ is equivalent to constructing a local exponential curve-fit to the arc 
length between the pre-assigned boundary grid points. The interpolation of 
the parameters into the interior of the computational domain simply extends 
the range of the curve-fit. The elliptic equation system (2.21) then merely 
provides a reliable, automatic means for translating the parameters into a local 
exponential curve-fit at each interior point that reflects the boundary value 
distribution, and that has the properties of regularity and monotonicity 
required of non-singular coordinate transformations. The resulting grid has 
the further desirable property that the two families of grid lines are locally 
orthogonal at the boundaries of the physical flow region. 

As a final observation, we point out that the general method has the flexibility 
to allow one to control at will the angle of intersection between the two families 
of grid lines at boundaries. This can be accomplished as follows. In place of 
the orthogonality condition (2.28); we use the generalized condition 

r ? |cos© on £ = (2.32) 

where 0 denotes the local angle of intersection between the boundary curve 
£ = Cb ar) d The family of transverse coordinate curves n = const. A more 
convenient representation of this condition is 


V r c 


= r_ 
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on c = % 


(2.33) 


y nW<; = ( Vr y c z n )cot 0 


which satisfies Eq. (2.32) identically. Upon inserting the zero-curvature 
constraint (2.26) into Eq. (2.25), all ^-differentiated terms in the resulting 
equation can be eliminated with the aid of Eq. (2.33). This yields the 
following equation for the parameter <p. 


= . 2 ( sin0 V) _ (yy ot0)y nn +(z n% cot0)z nn 


v +z 


This last equation can be used in the same fashion as Eq. (2.30) to compute 
numerically the parameter <p in terms of n-derivatives of the pre-assigned 
boundary values y,z and of any pre-assigned distribution of 0 as a function of 
position along the boundary curve £ = e;^. 
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Section 3 

EVALUATION OF NOZZLE THRUST AND DISCHARGE COEFFICIENT 


3.1 THRUST 

There are two methods that can be used to compute the nozzle performance 
characteristics from the converged steady-state flowfield solution. For 
example, the net thrust may be computed by integrating the sum of the axial 
components of pressure force and viscous shear stress over the surfaces of the 
nozzle wall. Alternatively, the thrust may be computed from a global momentum 
balance. We shall employ the latter method because it involves only the in- 
tegration of the Cartesian components of the momentum flux vector over the 
peripheral faces of the computational space; whereas the former method requires 
that computed velocity field be differentiated numerically in order to deter- 
mine the wall shear stress. 

The fact that a global momentum balance can be used to evaluate force components 
such as thrust is a formal consequence of the global conservation properties of 
the system of partial differential equations that govern the flow, namely, the 
Navier-Stokes equations (2.1). These global conservation properties can be 
derived simply by taking the volume integral of Eq. (2.1) over the entire com- 
putational space £, n » £• As a concrete example, let us consider the internal 
flow in an isolated three-dimensional nozzle with the right-handed curvilinear 
coordinate system £, n, z, defined so that the surfaces £ = £ Q anc * ^ = ^max 
represent the inflow and outflow boundaries, respect ively ; the surfaces n = n 0 

and n = n represent the left and right sidewalls, respectively ; and the 
max 

surfaces z, - z, Q and z, - represent the upper and lower walls, respec- 

tively. The volume integral of Eq. (2.1) over the computational space 

£ < £ < £ , n <n<n > c < £ < ? then may be written as 

^o — s — s max 7 o — - max o - - max J 
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d£ dn dc + 



d? dn dc 


+ 


+ 



dn dc d c 


dc dc dn 


0 


(3.1) 


The first term in Eq. (3.1) vanishes at steady state. The volume integral in 
each of the remaining terms reduces to the difference of two surface integrals 
because the argument of the volume integral is a perfect differential with re- 
spect to one of the three spatial coordinates. The final result may be written 
in the following form 



Re 


-1 


0 n max > c)] dc dc 


- J [gU> n 0 


0 - Re 1 eU, n Q > c)] dc dc 


f h(5. n. c max )-Re _1 aU, n, c max )] dC dn 


(3.2) 


/' 


[h(c> n, C Q )-Re ^ w (?> n, C Q )] dc dn 


-/ f(C 0 » n» ?) d n d ? / 


c) dn dc ~ n, c) dn dc 


max 


Equation (3.2) is a formal expression of the steady state global conservation 
principles for the entire computational space. If we examine the first scalar 
component of this vector equation, each of the terms on the left side of the 
equation represents the net flow of mass through one of the four nozzle walls. 
Similarly, in the fifth scalar component of (3.2), each term on the left repre- 
sents the net energy flux through one of the nozzle walls. In the middle three 
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components of Eq. (3.2), each term on the left represents the x, y, or z 
component of the total force acting on a wall. This includes the viscous 
shear force, which is represented by the term that is inversely proportional 
to the Reynolds number Re in the argument of the surface integral. 

The terms collected on the R.H.S. of Eq. (3.2) have a corresponding physical 
interpretation as the net flux of mass, momentum, and energy through the 
nozzle inlet and exit planes. Thus if we define a "generalized force vector" 

F whose five components represent, respectively, the net mass flux through 
the nozzle walls, the three Cartesian components of the net force acting on 
the nozzle walls, and the net energy flux through the nozzle walls; then F 
can be computed from the equation 

F =/( V ?) d n dc - n ’ d dr) d? (3.3) 

which is simply a restatement of Eq. (3.2). 


In the general case where both the interior flow and the external flow about 
a bilaterally symmetric nozzle are computed numerically, the nozzle structure 
is embedded entirely within the computational space. The curvilinear coordi- 
nate system then is defined such that the surface s; = represents the in- 
flow boundary and £ = £ max represents the outflow boundary, which is situated 
downstream of the nozzle structure; n = n Q and q = represent flowfield 
symmetry planes; and the surfaces n = n ^ and t, - w represent outer 
freestream boundaries. The steady-state volume integral of Eq. (2.1) over the 
computational space then yields the following equation for the generalized 
force vector F acting on the nozzle structure 


=/f( V d dn d ^ -/ 


d dn ds -/ fU_, n, 0 dn d£ 


max 


-fat’ "max’ d? d5 -fa ’ n> Snax^ d ^ dl1 


(3.4) 
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The viscous terms do not appear in the last two integrals on the right of this 
equation because those terms vanish in the inviscid freestream (Ref. 1, 

Section 2.4.4) . 

The surface integrals in Eqs. (3.3) and (3.4) can be evaluated numerically from 
the computed flow variables by sequential application of one-dimensional 
trapezoidal integration formulas for the two coordinate directions along the 
surface. The trapezoidal formula is appropriate because it is consistent with 
the second order spatial accuracy of the numerical algorithm that is used to 
compute the flowfield. 

We shall usually define the Cartesian coordinate system so that the x axis 
is oriented in the general streamwise direction. The second component of the 
generalized force vector F then represents the net thrust. The first and 
fifth components of the generalized force vector represent the net mass and 
energy fluxes at the walls of the nozzle. The net wall mass flux should be 
identically zero if the walls are impermeable. Similarly, for adiabatic wall 
boundary conditions , the net energy flux through the walls should also vanish 
identically. The extent to which these components of F differ from zero 
when evaluated numerically then provides a measure of the global accuracy of 
the flowfield computation. 

We observe that the described method of computing the generalized force vector 
from surface integrals over the boundaries of the computational space is valid 
only if the algorithm that is employed to compute the flowfield does not compro- 
mise the global conservation principles that are satisfied by the partial 
differential equations (2.1). That is, the difference equations derived from 
(2.1) using the algorithm must obey the same global conservation principles, 
or Eqs. (3.3) and (3.4) will be invalid. Great care has been taken in the 
numerical computation of boundary conditions and in the formulation of artifi- 
cial smoothing and dissipative operators to ensure that the composite numerical 
algorithm does indeed possess the same global conservation properties enjoyed 
by the original partial differential equations. 
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3.2 DISCHARGE COEFFICIENT 


The nozzle discharge coefficient C is defined as the ratio of the total mass 

w 

flow rate through the nozzle to the flow rate that would exist if the flow in 
the nozzle were isentropic [Ref. 4, p. 99]. According to one-dimensional 
isentropic nozzle flow theory, the dimensionless isentropic flow rate per unit 
area at the throat of a nozzle under flow conditions is [Ref. 4, p. 85] 



where y is the specific heat ratio of the gas and the dimensionless density 
and velocity are referred to the stagnation chamber density and sound speed, 
respecti vely. The discharge coefficient is then given by the equation 



where Wy denotes the actual mass flow rate at the nozzle throat and Ay is 
the cross-sectional area of the nozzle at the throat. The latter two quanti- 
ties may be computed from the final steady-state flowfield solution as follows. 


We assume that the curvilinear coordinate transformation has the form £ = £(x) , 
so that the surfaces £ = const, represent cross-sectional planes. Let 
£ = £y denote the throat location in the curvilinear coordinate system; the 
surfaces n = n Q and e = represent flowfield symmetry planes; and the 
surfaces n = n w and c = c w represent complementary portions of the 
interior surface of the nozzle wall. The first component f-j of the flux 
vector f can be interpreted physically as themassflux in the 5 coordinate 
direction. The total mass flow rate at the throat then is simply 


W 


T 



n» c) dri dc, 


(3.6) 
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An analogous expression for the throat cross-sectional area can be obtained by 
noting that the normal vector to any surface £ = const, has the Cartesian 
representation 

n = V£ (3.7) 


whence the element of area in such a surface is [Ref. 5] 
dA = j| vel = I /?* + T 2 + T\ 

Thus, the throat cross-sectional area is simply 

p w ( — 

/ \jl\ + l 2 + 5* dn dc 

•'o 



(3.8) 


(3.9) 


where the metrics £ , £ , £ are evaluated at the throat location £ = £ T 

The surface integrals that appear in Eqs. (3.6) and (3.9) can be evaluated 
numerically in the same fashion as outlined in the preceding subsection for 
the surface integrals in the equations for the generalized force vector F. 
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Section 4 

FINAL FORMULATION OF TURBULENCE MODELS 

The viscous terms on the R.H.S. of Eq. (2.1) involve the Prandtl number and a 
non-dimensional viscosity coefficient p > which is referred to the dimensional 
viscosity at some reference state at which the Reynolds number Re is defined. 
For laminar flow, these are molecular transport properties, which we shall de- 
note by the subscript e. For air at moderate temperatures, the Prandtl 
number is approximately constant 

Pr e = 0.72 (4.1) 

and the variation of viscosity with temperature may be approximated by a power 
law 


u e - T m (4.2a) 

or may be computed from the Sutherland law 

t 3 /2 

p e = T + 110.3/T r (°K) (4.2b) 

where the dimensionless temperature T is referred to the dimensional reference 
temperature T r at which the Reynolds number Re is defined. 

For turbulent flow, the viscosity is taken as the sum of the molecular value 
and a turbulent eddy viscosity 

p = u e + pf (4.3a) 

where the dimensionless eddy viscosity is normalized by the reference value of 
u e that is used in evaluating the Reynolds number. The thermal conductivi ty. 
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which is proportional to the ratio between yi and P r, is also taken as the 
sum of laminar and turbulent contributions 


v_ 

Pr 




(4.3b) 


where and Pr^ are obtained from some sort of turbulence model. Such 
models generally employ a constant value for the turbulent Prandtl number 


Pr t = 0.9 (4.4) 

whereas the eddy viscosity y 0 is strongly dependent on the character of the 
flow; e.g., boundary layer, shear layer, wake, or jet. The models that are 
used here for nozzle flows are presented below. The models are discussed 
first in the context of two-dimensional and axisymmetric flows, and are later 
generalized to more complicated three-dimensional flows. 


4.1 TWO-DIMENSIONAL AND AXISYMMETRIC FLOW 

In the present application, we require a turbulence model that is valid in the 
nozzle wall boundary layers; in the wake region behind a nozzle wall, side- 
plate, or wedge-plug; in the near-field mixing layer between the external flow 
and the nozzle exhaust stream, and in the far-field fully-developed jet region. 
Standard engineering turbulence models for the eddy viscosity are restricted to 
one or another of the described sub-regions of the flowfield, and must be 
patched together to provide a composite model. A general discussion of 
engineering models that apply in the various sub-regions has been given in 
Section 6 of Reference 1. These standard models were designed originally for 
use in analytical or simple numerical solutions for flows where the turbulent 
region is essentially a two-dimensional thin layer adjoining a region of 
spatially uniform, inviscid flow. In this type of solution, the boundaries 
of the turbulent region are relatively well-defined (such as the wall and the 
outer edge of a wall boundary layer) and determine the length scale in terms 
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of which the turbulence model is formulated. In complicated numerical solu- 
tions such as the nozzle flows of interest here, this need to locate the 
physical edge of a boundary layer, shear layer, or wake poses considerable 
difficulty; both because there generally are substantial flow gradients eyen 
in inviscid regions, and because spatial and temporal oscillations often exist 
in the computed flow variables at grid points. The difficulty, already severe 
in two-dimensional flows, can become extreme in three-dimensional flows where 
the geometry of a shear layer, for example, may be so complicated that it is 
not clear how one ought to proceed in order to find the boundaries of the 
layer and use that information to define a local length scale. Furthermore, 
computational accuracy can be highly uncertain because the engineering turbu- 
lence models generally are quite sensitive to the numerical value of the local 
length scale. 

Baldwin and Lomax [Ref. 6] recently have presented a turbulence model for two- 
or three-dimensional wall boundary layers and wakes that does not require 
finding the boundaries of the turbulent region. This model is based directly 
on the engineering turbulence models described in Section 6 of Reference 1, but 
uses the spatial distribution of vorticity to determine the length scale in 
terms of which the eddy viscosity is computed. We shall employ a modified 
version of this model for the nozzle wall boundary layers and for the near- 
wake region downstream of the trailing edge of a nozzle wall, side-plate, or 
wedge-plug. For mixing layers and for the fully-developed jet region, we have 
developed a simple Prandtl mixing length type of model in which the turbulent 
length scale is defined in terms of the vorticity distribution, rather than 
in terms of the physical width of the mixing layer or jet. To facilitate the 
description of the models, we first define the character of the curvilinear 
coordinate system that is used for two-dimensional or axi symmetric flow. 

We orient the x-axis of the Cartesian base coordinate system in the general 
streamwise direction. For two-dimensional flow in the x-z plane, the flow is 
invariant with respect to y. The right-handed boundary-conforming curvilinear 
coordinate system is defined such that g = £(x), n = n(yK c = c(x,z). For 
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ax i symmetric flow, the curvilinear coordinates are defined as in Eqs. (2.12), 
so that n is the azimuthal coordinate. The turbulence models to be used in 
the various sub-regions of the flowfield are presented below in terms of the 
described coordinate system. For convenience, the models will be presented in 
terms of a turbulent kinematic viscosity v t that is scaled by the reference 
Reynolds number. The dimensionless eddy viscosity coefficient that enters into 
Eq. (4.3) is obtained from v t by the following equation 

u t = P Re v t (4.5) 


4.1.1 Wall Boundary Layers 


The Baldwin-Lomax turbulence model [Ref. 6] is a two-layer model in which 
is given by 


(v t ) 


inner 


K) 


outer 


4 < 4„ 
— c 


4 > 4 


V t 


(4.6) 


where 4 denotes the normal distance from the wall and 4„ is the least value 

c 

of 4 at which the inner and outer viscosities are equal. The viscosity for 
the inner region is defined by 



= £ 2 \Z\ 

(4.7a) 


£ = fe4 [1 - exp (-4 + /A + )] 

(4.7b) 


4 + = — \l Re p t 

u ew V h w w 

(4.7c) 

where 

S = vx V 

(4.8) 
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is the vorticity, the subscript w denotes conditions at the wall, the wall 
shear stress is given by 


T w y ew 


and the constants are 


k = 0.4 

A + = 26 


The corresponding formula for the outer region is 


(v t ) 


outer C 1 F F k^ 


C- ( = 0.0269 


(V-6 f 

I mm 

I* f 

V m 


K < 1 


K > 1 


K = U/2f 

m 

where U is the maximum velocity in the profile 
U - |V(i)l 

F^U) is the Klebanoff intermittency factor 
F k U) = [1 + C 2 (C 3 A/^) 6 ]" 1 

C 2 = 5.5 C 3 = 0.3 


(4.9) 


(4.10) 


(4.11a) 

(4.11b) 

(4.11c) 
(4. lid) 

(4.12) 

(4. l 3a) 
(4.13b) 
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and where the quantities 4 m and f are defined at the maximum point of the 
function 


fU) = 4|S| [1 - exp (-//A + )J (4.14) 

In a numerical solution where both 4 and f are known only at discrete grid 
points, Baldwin and Lomax [Ref. 6 ] recommend that the true maximum be found by 
employing a three-point quadratic fit to the function f(4) in the neighborhood 
of the apparent maximum that occurs on the grid. 

This model has been used successfully in the computation of two- and three- 
dimensional flows with either attached or separated boundary layers [Ref. 6 ]. 

To apply the model to two-dimensional and axisymmetric flows in the curvilinear 
coordinate system n, ? described earlier, we make several approximations 
that sharply reduce the computational labor. First, we assume that the £ 
coordinate is approximately orthogonal to the wall, so that the normal distance 
4 can be approximated as the arc length along ? coordinate curves. If r 
denotes the radius vector in the Cartesian coordinate system, then the vector 
that is locally tangent to a 5 coordinate curve is 

- (x c , y^, z ? ) (4.15) 

The elemental arc length along the curve is then 

^ = |rj = ^ + y* + z* (4.16) 

and the distance along the ? coordinate curve can be obtained by integration 
of Eq. (4.16) from the wall outward. 

The second approximation that we shall use relates to the computation of the 
vorticity vector in Eq. (4.8). The latter can be expanded in terms of the 
curvilinear coordinate system by means of the chain rule to obtain the equiva- 
lent expression 
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0 ) 


V x V 


( 4 . 17 ) 


“)■ 

= V£ x \L+ Vr) x V + V£ x V 
£ n C 

For a two-dimensional flow in the surface n = const., one can show easily that 
the dominant contribution to the magnitude of the vorticity vector in Eq. (4.17) 
comes from the component that is normal to the surface n = const. This state- 
ment also holds for axisymmetric flows in which n represents the azimuthal 
coordinate. The unit normal vector to a surface n = const, is 


n = vn/ 1 Vti I 


( 4 . 18 ) 


and the vorticity component in this direction is the inner product 


n 




(4.19) 


By invoking the well-known vector identity governing successive inner and 
cross products, (4.17) and (4.19) can be combined to yield the result 


~y ~y 

n * w 


[(vn x vs)'V ? + 


(vn x v?) • V ]/ | vn I 


(4.20) 


The first term in brackets in this last equation involves quantities that we 
already have neglected in deriving the parabolized Navier-Stokes equations 
(2.1) [Ref. 1, Section 2.3], and we shall neglect this term here as well. In 
fact, one can see from Eq. (4.20) itself that the first term involves velocity 
derivatives with respect to the streamwise coordinate £ ; i.e., in the di- 

rection along the wall. Such derivatives are always small in a boundary layer 
compared to derivatives V in the direction away from the wall, and can 
safely be neglected. With this approximation, the dominant part of the vorti- 
city magnitude can be written in the very simple form 


• V 

JL X 

J|vn| 


X U c v c + w t ! 

+ + "I 


( 4 . 21 ) 


where we have made use of the identity 
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vn x vc 


(4.22) 


= J 


-1 


which is a property of the coordinate transformation x,y,z n, c. 


The numerical implementation of the above turbulence model is accomplished as 
follows. The spatial derivatives in Eqs. (4.16) and (4.21) are approximated 
by central difference operators at interior grid point and by the appropriate 
forward or backward difference operator at boundary grid points. The arc 
length a is computed from 6^ by using the trapezoidal quadrature formula. 
Finally, the true maximum point of the function f in Eq. (4.14) is found as 
a function of the c coordinate; i.e., by operating in the computational co- 
ordinate system (s, n, c) in which the grid spacing is always uniform, rather 
than in physical space where the grid spacing is nonuniform with respect to 
the arc length 6. The c coordinate will always be stretched with respect to 
the physical coordinate .6 in order to resolve the flowfield gradients that 
exist in the boundary layer, and a local quadratic fit to the function f(s) 
is performed more easily and more accurately than a corresponding fit to the 
function fU). Once the maximum point ? is found, a similar quadratic fit 
to the function a(c) is used to determine the corresponding value a = i(? m ). 


4.1.2 Wakes 


Baldwin and Lomax [Ref. 6] state that the formulation in Eqs. (4.11) - (4.14) 
for the outer region of a boundary layer also can be used in wakes if the 
bracketed exponential factor is omitted from Eq. (4.14), and U is redefined 
as the difference between the maximum and minimum velocities in the wake region 


U 



l v U 


(4.23) 


Although it is not so stated in Reference 6, the transverse coordinate 6 
presumably is reckoned from the point of minimum velocity. For an asymmetric 
wake, this would imply that the regions on either side of the velocity minimum 
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are to be treated separately. If this were done, it follows from Eqs. (4.11) 
and (4.14) that the turbulent viscosity would be discontinuous at the point of 
minimum velocity. Furthermore, Baldwin (in a private communication) has stated 
that computational instabilities were encountered in attempting to use the outer 
formulation alone in the near wake of an airfoil because the grid is most re- 
fined near the velocity minimum where the value of is greatest. The insta- 

bility can be avoided (according to Baldwin) by arbitrarily using the same two- 
layer formulation as in the boundary layer, with the bracketed exponential factor 
omitted from Eq. (4.7b). Because the resulting inner formulation is inconsistent 
with other wake turbulence models, we shall employ a different model for the 
inner region of the wake, i.e., near the velocity minimum. This model employs 
the Prandtl mixing length formula (4.7a) but the mixing length l is defined as 




C U/IJ 


max 


(4.24) 


where U is given by Eq. (4.23), C is a constant and 


max 


max 
-6 >0 


|SU) 


is the maximum vorticity magnitude in the section of the wake under consideration, 
since the two sections on either side of the velocity minimum are treated 
separately. Note that v t in Eq. (4.7a) remains virtually continuous at the 
border between wake sections even though the latter are treated separately, be- 
cause |S| as computed from Eq. (4.21) essentially vanishes at the velocity 
minimum. 


For a wake, the constant C in Eq. (4.24) has the value 


C 


wake 


0.255 


This value was obtained as follows. 


(4.25) 


In any turbulence model such as that of Eqs. (4.7a) and (4.24), the constants 
must be evaluated from experimental data for the specific type of flow under 
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consideration. This is done for wake turbulence models by requiring that the 
wake spreading rate predicted by the model match that measured experimentally. 

An analytical wake solution based on a mixing length model is given in Schlichting 
[Ref. 7, p. 600]. In this solution, the mixing length is defined differently 
from Eq. (4.24), and is evaluated so as to match the measured wake spreading 
rate behind a circular cylinder transverse to the flow. We have used the wake 
velocity profile from the analytical solution to deduce the value of C in 
Eq. (4.25) that matches the experimental spreading rate when the mixing length 
is defined by Eq. (4.24). 

Schlichting [Ref. 7, p. 603] also gives an analytical wake solution based on 
a different eddy viscosity model in which is given by 

v t = k b U (4.26) 

where b is the wake width and < is a constant that is determined to match 
the measured wake spreading rate. This model, unlike the mixing length model, 
assumes to be constant throughout the wake. Yet, the velocity profile 
from this analytical solution can be used to deduce a value for the constant 
C in the mixing length model that is in very close agreement with the value 
in Eq. (4.25). The value C wa(<e = 0.259 is obtained by simply requiring that 
the mixing length model yield the same value for at the maximum vorticity 
point in the wake as does the model (4.26) from which the analytical velocity 
profile is derived. This means that we can calibrate the Prandtl mixing length 
model for a given type of flow by using the analytical solution from the simpler 
model (4.26) once the latter has been calibrated to agree with experiments for 
that type of flow. We shall make use of this calibration technique below in 
modeling for mixing layers and jets. 


4.1.3 Mixing Layers and Fully-Developed Jet Region 

As discussed in Section 6 of Reference 1, we had intended to employ the constant 
eddy diffusivity model (4.26) in .mixing layers and in the fully-developed region 
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of the exhaust jet. This type of model in which is a constant is diffi- 
cult to link to other models for adjoining flow regions without having large 
discontinuities in the spatial distribution of v^.. For this reason, we 
favor the Prandtl mixing length model given in Eq. (4.7a). There, v-j. 

is proportional to the local vorticity, which tends to be small both in inviscid 
regions and near the border between turbulent sub-regions that require different 
turbulence models. We shall apply a one-layer model of the type (4.7a) for a 
mixing layer or jet. As in the inner region of a wake, we determine the mixing 
length scale from the velocity extrema and the maximum vorticity (Eq. (4.24)), 
but the constant C must be calibrated separately for each type of flow. 


Mixing Layers 

Schlichting [Ref. 7, p. 598] gives an analytical solution for the two-dimensional 
mixing layer based on the constant eddy viscosity model (4.26), with the con- 
stant k selected to match experimental data on the width of the turbulent 
region. We have calibrated the mixing length model by the technique described 
in Section 4.1.2; that is, by requiring that the mixing length model yield the 
same value for at the maximum vorticity point in the layer as does the 

constant eddy viscosity model. The resulting value of the constant C in 
Eq. (4.24) is 

C . = 0.136 (4.27) 

mix 


Fully-Developed Jet 

The region downstream of the point where the inner edge of the mixing layer be- 
tween the external flow and the nozzle exhaust stream penetrates to the flow 
centerline is known as the fully-developed jet. Since we are interested pri- 
marily in three-dimensional flows, we assume that this region far downstream of 
the nozzle exit is essentially similar to that for a round (axi symmetric) jet. 
Accordingly, we calibrate the mixing length model from the analytical solution 
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for the round jet based on the constant eddy viscosity model [Ref. 7, p. 607] s 
by again requiring that the mixing length model yield the same value for 
at the point of maximum vorticity. This gives the following value for the 
constant in Eq. (4.24) 

C jet = 0.129 (4.28) 

The latter is very close to the value in Eq. (4.27) for a mixing layer. Since 
the turbulence models are only approximate, in programming the models for numeri- 
cal solution, we simply have used the same value 

C = 0.13 (4.29) 

for both mixing layers and jets. 


4.2 GENERAL THREE-DIMENSIONAL FLOWS 

For axisymmetric three-dimensional flows, the turbulence model equations given 
in the preceding section are obtained by viewing the flow as essentially two- 
dimensional in an azimuthal coordinate surface n=const. We denote the turbulent 
viscosity so obtained as v> , since Eq. (4.21) for the vorticity magnitude in- 
volves V in which derivatives of the velocity are taken with respect to the t, 
coordinate direction within a surface n = const. For general three-dimensional 
flows governed by the parabolized Navier-Stokes equations (2.1) in which the 
viscous terms associated with the streamwise coordinate £ are neglected, there 
are two principal cross-stream coordinate directions n, c- In this case, we use 

(t) 

the same quasi -two-dimensional approach to obtain ' based on regarding the 
flow as two-dimensional within a coordinate surface n = const., and apply a 
similar quasi-two-dimensional approach to compute a second value based 

on regarding the flow as two-dimensional in a surface c = const. The equations 
for v[ n) are obtained from Eqs. (4.6) - (4.29)by the substitution (n, c)-*U, n). 
To obtain a single composite value for at each point of flowfield, we arbi- 
trarily combine the quasi-two-dimensional values by using the root-mean-square 
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X 

2 


( 4 . 30 ) 


V 


t 




In mixing layers and jet regions, both and are obtained from 

the same formulation. However, in wall boundary layers and near wake regions, 
a certain arbitrariness exists in the quasi-two-dimensional approach. Con- 
sider, for example, a wall boundary layer. The boundary-conforming character 
of the curvilinear coordinate transformation is such that a wall is represented 
as either a surface n = const,, a surface z, = const, or as a composite of 
intersecting surfaces of each type. Over any portion of a wall that is repre- 
sented as a surface c = const., the nature of the transformation is such that 
coordinate lines along which c alone varies are nearly orthogonal to the 
wall in the boundary layer, and the boundary layer formulation given in Eqs, 
(4,16) - (4.21) applies for computing . Within each wall-like coordinate 

surface c = const., the boundary layer formulation (4.6) - (4.21) with the 
substitution (n,c) ->■ (?, n) also is used to compute in boundary layer 

regions where this surface ? = const, intersects a wall that is represented 
as a member of the other family of coordinate surfaces n = const. However, 
a different model is necessary when this is not the case. An example is an 
axi symmetric flow where n is the azimuthal coordinate and the wall is a sur- 
face s = const, [cf. Eq. (2.12)]. The boundary layer formulation then is 
appropriate for but is inappropriate for vj. 1 ^ since the c = const, 

surfaces are wall-like coordinate surfaces that do not intersect the wall. In 

such cases, we arbitrarily use the mixing layer formulation to compute vi n ^ . 

/ \ t 

This always yields a value of vl r]) that is dominated by the value of 

(r ) z 

v> 'obtained from the boundary layer formulation, and ensures that the composite 

z (c) 

value in Eq. (4.30) will recover the correct result, namely, vj.' , within 

the boundary layer. The same approach is used in two-dimensional wake regions 

in that the quasi-two-dimensional value of v t for the coordinate surfaces 

transverse to the wake is obtained from the mixing layer formulation. 
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SECTION 5 


NUMERICAL EXPERIMENTS 


A number of numerical experiments have been performed for two-dimensional internal 
and external flows to test various aspects of the implicit numerical method. These 
experiments are described in the following subsections. Section 5.1 demonstrates 
that the numerical method is capable of computing internal flows with boundary 
layer separation. Section 5.2 gives the results of computations that have been 
performed to investigate the effects of implicit boundary conditions and of time 
stepsize on the rate of convergence to steady state and on the final steady state 
solution. Section 5.3 gives the results of parametric calculations to investigate 
the sensitivity of the solution to the magnitude of the explicit smoothing co- 
efficient. Section 5 demonstrates that the artificial implicit dissipative terms 
described in Section 2.4 have a favorable effect on the numerical stability and 
on the rate of convergence to steady state. The major conclusion drawn from these 
experiments are summarized in Section 5.5. 

5.1 COMPUTATION OF SEPARATED FLOW 

To demonstrate that the numerical method is capable of computing internal flows 
with boundary layer separation, we have performed an internal flow computation for 
a two-dimensional converging-diverging nozzle whose wall is composed of two cosine- 
shaped segments. The nozzle shape and computational grid are shown in Fig. 5-1. 

The computation is for laminar flow at a Reynolds number based on stagnation 
conditions of 10 s , a Prandtl number of unity, and viscosity proportional to temp- 
erature. The grid is 15 x 15 in the x and z directions, stretched exponentially 
in the z direction to resolve the nozzle wall boundary layer. That is, the 
curvilinear coordinate transformation is defined by the equations 

z/ z maxU) = [exp(ac/c max )-l]/[exp(a)-l] 

C ~ ( &“ 1 ) AC , AC — 1 , 1,2,..., 1 

where a is a constant stretching coefficient, and z max (x) is the wall shape. 
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Figure 5-1 Computational Grid for Cosine Nozzle 


Crude a priori estimates of the boundary layer thickness indicated the latter to 
be approximately one percent of the local nozzle height over most of the length 
of the nozzle wall. To resolve this thin region, the stretching parameter a was 
selected to produce a grid spacing at the wall that is a factor of 5000 smaller 
than the grid spacing at the centerline. 

The computation has been performed using the implicit inflow, outflow, symmetry, 
and adiabatic wall boundary condition computation schemes described in Sections 
2.4 and 4.2 of Reference 1 and in Section 2.3 herein. 

The initial conditions are as follows. The nozzle is assumed to be choked at the 
throat. The pressure, density, temperature, and streamwise velocity (averaged 
over the cross section) are computed from one-dimensional inviscid isentropic flow 
theory for the nozzle area variation. These inviscid core conditions are applied 
over the lower half of the grid 0<c<£max/2> assuming that the streamlines are 
parallel to the grid lines £ = constant. The velocity components at the nozzle 
wall £ = £ max are set to zero to satisfy the no-slip viscous wall boundary conditions. 
The velocity components at the remaining grid points are linearly interpolated in 
C between C max /2 and £ max - T ^e Crocco relation is used to compute temperature 
from velocity in this region. Density follows from the equation of state by taking 
the pressure as uniform over the nozzle cross section at each x station. 

Convergence to steady state was obtained in 300 time steps using a constant step- 
size At = 0.05, which corresponds to a Courant number Co(0 = 530. The computed 
longitudinal pressure distributions along the wall and centerline are shown in 
Fig. 5-2, together with the pressure distribution predicted by one dimensional 
isentropic flow theory (Ref. 4) for the nozzle area variation. Pressures shown 
are referred to stagnation pressure. As one might expect, the computation shows 
a substantial recompression at the wall near the exit plane that is induced by the 
locally concave wall shape. The inviscid core flow is supersonic in this region, 
and the adverse pressure gradient causes boundary layer separation between x * 1.4 
and the exit plane. The longitudinal velocity profile in the separated region is 
smooth, as shown in Fig. 5-3, and displays the classic flow reversal in the near- 
wall region that one expects from boundary layer theory. In the figure, velocity 
is referred to the stagnation sound speed. 

The vertical profiles of pressure and velocity at the geometric throat, x = 0, are 
given in Figs. 5-4 and 5-5, and show substantial nonuniformities due to the two 
dimensionality of the flow. 
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Figure 5-2 Longitudinal Distributions of Pressure, Cosine Nozzle, Re 
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Figure 5-3a Velocity Profiles Across Cosine Nozzle in Separated Region, 

Re = 10 s . 


47 



NOZZLE WALL 



Z/ZMflX 


CYCLE 300 TIME-15.000 K- 1 J- 8 



P 


Figure 5-4 Pressure Profile Across Cosine Nozzle at Throat, x = 0. 
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Figure 5-5 Velocity Profiles Across Cosine Nozzle at Throat, x - 0. 
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5.2 EFFECT OF IMPLICIT BOUNDARY CONDITIONS AND TIME STEPSIZE ON THE SOLUTION 


5.2.1 Implicit Boundary Conditions 

A set of numerical experiments has been performed to explore whether the use of 
implicit boundary conditions has any significant effect on the rate of convergence 
to steady state or on the final steady state solution. These tests were made for 
two-dimentional laminar flow over an adiabatic flat plate at a freestream Mach Number 
= 0.1, Reynolds Number Re^ = 10 5 , Prandtl Number Pr = 1, and viscosity propor- 
tional to temperature. The grid and initial conditions are as described in Ref. 1, 
Section 5. In one case, implicit wall, freestream, and outflow boundary conditions 
were employed as described in Ref. 1, Sections 2.4 and 4.2. In the other case, 
time-lagged explicit boundary conditions were applied as follows. During each time 
step, the flow variables q n+1 = q n +Aq are computed implicitly at interior grid points 
by assuming that Aq = 0 at boundary grid points. The flow variables at boundary 
points are then updated at the end of the step by extrapolation of the solution 
q n+1 f rom interior points. At wall boundary points, the velocity components (u,v,w) n+ * 
are set to zero and the wall pressure is computed from the normal momentum equation. 
For a flat plate, the latter equation implies that the pressure gradient in the 
wall-normal direction vanishes, hence the pressure at each wall point is set equal 
to that at the adjacent interior point. To satisfy the adiabatic wall conditions 
the gas temperature at the wall is extrapolated quadratical ly with zero gradient 
using the temperature at the two nearest interior points along the wall-normal grid 
line. At the downstream outflow boundary, q n+1 is extrapolated linearly from the 
two nearest interior points along each grid line transverse to the outflow boundary. 

At the lateral outer boundary, the freestream values of pressure, temperature, and 
streamwise velocity component u are imposed, and the transverse momentum flux 
components pv and pw are extrapolated linearly from interior points. 

Numerical solutions for the two cases, implicit vs. time-lagged boundary conditions, 
showed no significant difference in either convergence rate or in the final steady- 
state solution as obtained with a constant time stepsize At = 0.01. The latter 
corresponds to a Courant number Co = 40 based on the minimum mesh spacing in the 
direction normal to the plate. Both cases converged within 300 steps and yielded 
steady-state drag coefficients that differed by only 0.07%. However, attempts 
to increase the time stepsize revealed that the solution with implicit boundary 
conditions remains numerically stable with time stepsizes significantly greater 
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than those attainable with time-lagged boundary conditions. 

With implicit boundary conditions, the computation remains stable and converges 
well for a time step as large as At = 0.1 (Co = 400). With time-lagged boundary 

conditions, the computation is unstable for At - 0.05 (Co = 100). We have not 
attempted to determine the precise stability boundaries for the two cases, because 
our experience has shown that the stability boundary depends upon the flow 
conditions. With Moo = 3, for example, the flat plate flow computation is unstable 
at At = 0.05 (Co = 300), even with implicit boundary conditions. 

As shown in the next section, the solution usually converges to steady state more 
rapidly, i.e., in fewer steps, the larger the stepsize At. Because the use of 
implicit boundary conditions permits a larger stepsize, we conclude that implicit 
boundary conditions significantly improve overall computational efficiency by 
allowing convergence to be attained in fewer steps than would be required with 
time-lagged boundary conditions. 

In addition to the described experiments comparing implicit and time-lagged boundary 
conditions, we have performed some test computations to evaluate the sensitivity 
of the solution to the placement of the lateral outer boundary at which freestream 
pressure, temperature, and streamwise velocity component, u, are imposed as boundary 
conditions, but the transverse momentum components are computed implicitly from 
the transverse momentum equations. The results of these tests have been reported 
in detail in Ref. 9, and we merely summarize them here. 

For supersonic flow, ^>1, the solution was found to be insensitive to the location 
of the outer boundary as long as the latter was positioned far enough above the - 
plate to enable the shock wave generated by the viscous interaction to be captured 
in the mesh. 

For subsonic flow, Moo<l, the solution again was insensitive to the placement of 
the outer boundary. The overall solution and the computed drag coefficient differed 
by less than 0.1% when the outer boundary was shifted inward from ten boundary layer 
thicknesses above the plate to as little as two boundary layer thicknesses above' 
the plate. The convergence histories were also virtually identical for the two cases. 
Similar computations in which all freestream conditions were imposed at the lateral 
outer boundary displayed a markedly different behavior. As much as a factor of 
three more time steps were required to attain convergence to steady state, and the 
final solution was quite sensitive to the placement of the outer boundary. The 
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computed drag coefficient was found to vary with boundary position, and the 
velocity profiles within the boundary layer showed a serious distortion that did 
not arise when the transverse momentum components were computed implicitly at 
the outer boundary. 

5.2.2 Time Steps ize 

As indicated in Ref. 1, a linear stability analysis shows that the implicit 
numerical algorithm is unconditionally stable for arbitrarily large values of the 
time step At. In practice, we have found that for a given type of flow, there is 
a limiting value of the time step beyond which the algorithm is unstable. The 
instability is attributed to nonlinear effects that are unaccounted for by linear 
stability theory. The evidence suggests that the instability may be associated 
with the streamwise coordinate direction ?(x) in which viscous effects are small 
and have been neglected in the governing equations (2.1). This hypothesis is 
based on the following observations. 

Let Co(?) denote the Courant number based on the inviscid flux vector f for the 
streamwise ? direction 

Co(?) = | A | max At/A? (5.1) 

A 

where A is an eigenvalue of the Jacobian matrix F in Eq. (2.7d) and where the 
maximum is taken over all five eigeuvalues and over all grid points. The Courant 
numbers Co(n) and CoU) associated with the transverse directions n»C are 

A A 

similarly defined in terms of the eigenvalues of the Jacobian matrices G and H. 

One can show that most explicit algorithms for the Navier-Stokes equations are 
subject to the stability criterion 

max{Co(?) , Co(n) , Co(?)> £ 1 (5.2) 

as long as the grid is locally fine enough to resolve the steep flow gradients 
that exist within viscous regions such as wall boundary layers and mixing layers. 
Note that the transverse Courant numbers are usually much greater than the stream- 
wise one because of the much smaller transverse physical grid spacing needed to 
resolve the steep transverse gradients that occur in viscous regions. 
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Co(n), Co(c)»Co(0 


(5.3) 


In contrast, the implicit algorithm remains stable for Co(n), Co(^)»l. However, 
with no artificial implicit dissipation (see Sections 2.5 and 5.4), we have found 
that the implicit algorithm tends to become unstable when the streamwise Courant 
number is of the order of unity or greater 

[co(C)J U n S tabl e ~ 0(1) (5.4) 

This still gives the algorithm a great advantage over explicit schemes. For 
relatively simple external and internal flows, one can obtain a crude estimate 
of the stable range of At from Eq. (5.4) when E, = £(x) by using the approximation 

Co(£) ~ c(1+M)At/Ax (5.5) 

where Ax is the physical grid spacing in the streamwise x direction and c and M are 
the dimensionless local sound speed and Mach number, respectively. For example, 
for external flows such as the flat plate boundary layer addressed earlier, the 
flow variables are referred to freestream conditions and Eq. (5.5) can be evaluated 
at the freestream boundary c = c 00 =l,M = M 00 using the smallest grid spacing 
(Ax) m ^ n « For internal nozzle flows where the flow variables are referred to 
stagnation condi tions, the streamwise grid is usually finest near the throat where 
M “ 1, c =-y/2/b+l) . 

External Flow Experiments 

We have performed a few numerical experiments to investigate the effect of the 
time stepsize on the rate of convergence for the = 0.1 adiabatic flat plate 
boundary layer problem described in Section 5.1.1. For a 15 x 15 grid, the 
stable range of At as estimated from Eq.'s (5.4) and (5.5) is 

(^stable ~ 0.065 

as noted in Section 5.1.1, the solution has been found to be stable for At as 

large as O.l|co(c) = 400 j and unstable for At = 0.s|co(c) = 2000 J. Table 5.1 

summarizes convergence data obtained from two computations with stepsizes that 
differed by a factor of five within the stable range. The first column in the 
table gives the stepsize, and the second column gives the time step number n at 

which the data in the remaining columns apply. The third column gives the 1_2 
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residual. R, which is defined as the volume integral over the computational space of 
the square of the set of steady-state terms in the Navier-Stokes equations (2.1) 
and the smoothing terms, normalized by the total volume 


R = (5v) 




d 5 d n d c (5.6) 


■! 


v = d^d n d^ 


The factor of one-fifth is introduced so that the residual represents an average 
over the five scalar components of the vector equation (2.1). The fourth column 
in the table lists the maximum relative change over the n'th step that occurred 
in any of the principal flow variables p, pu, and e at any point in the grid. 
During the later stages of convergence, this turns out to be the variable pu at 
one of the interior grid points adjacent to the wall. Note that this quantity 
strongly affects the local velocity gradient and hence the computed wall shear 
stress. This means that the computed skin friction drag represents a very 
sensitive indicator of whether the solution actually has attained convergence. 
The data in the last two columns of Table 5-1 are derived from this value 
| Aq/q | max - The fifth column contains an estimate of the instantaneous degree of 
unsteadiness in the solution 


IVjmax = |(Aq/AT)/q| max 


= i^Aq/q 


max 


(5.7) 


whereas in the last column, this quantity is scaled by the square of the time 
step. 


From Table 5-1, one can see that, after a given number of time steps, the solution 
for the larger time step is much nearer convergence. All three measures of 
convergence, the residual, the maximum relative Aq, and the degree of unsteadiness 
q T /q, are much smaller for At = 0.05 than for At = 0.01. The entries in the last 
column of the table are nearly independent of the stepsize for a fixed number of 
time steps n. This implies that, during the late stages of the calculation where 
the solution is near steady state, the remaining degree of unsteadiness | qx/q | max 
is inversely proportional to the square of the stepsize. We infer that the 
convergence rate is quadratic in the stepsize At, and that convergence is attained 
much more rapidly the larger the stepsize. 
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Convergence Data for Adiabatic Flat Plate Boundary Layer, 
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Internal Flow Experiments 

The time stepsize also has a strong effect on the convergence rate for internal 
flows. As an illustration, we present some results for a two-dimensional compu- 
tation of internal laminar flow in the vertical symmetry plane of a so-called 
"two-dimensional" converging-diverging nozzle that has been tested experimentally 
in the static test facility of the NASA Langley 16-foot Transonic Wind Tunnel. 

The nozzle configuration is described in detail in Section 6.1. We note here only 
that the actual nozzle has straight sidewalls and a rectangular cross section of 
constant width. The inlet section is of constant height. This is followed by a 
straight-walled converging section that is connected to a straight diverging section 
by a circular arc that forms the geometric throat region. 

The present two-dimensional flow test computation is for nozzle operating conditions 
corresponding to the nozzle design condition with a stagnation pressure of 1 atm. 
and a stagnation temperature of 295 Kelvin. The Reynolds number based on stagnation 
chamber conditions and throat half-height is 930,000. Although one would expect 
turbulent flow at this large a Reynolds number, the computation assumes laminar 
flow with a Sutherland viscosity law, a Prandtl number of 0.72, and adiabatic wall 
boundary conditions. The flow is computed in the upper half of the vertical plane 
of symmetry y = 0, z ^ 0, which is covered by a 23x15 grid in the x(j) and z(l) 
directions, respectively. The vertical (z) grid is stretched exponentially to 
resolve the nozzle wall boundary layer. The initial conditions used to start 
the computation were obtained in the same fashion as for the cosine nozzle test 
case discussed in Section 5.1. 

Two computations were performed, each starting from the same initial conditions. 

The first case was run for 500 steps with a constant stepsize Ax = 0.05, which 
corresponds to a Courant number of 500 based on the minimum grid spacing in the 
z direction. The solution at step 500 is only partially converged, as can be seen 
from Figs. 5-6 and 5-7, which show the first two components of the "generalized 

~y 

force vector" f (see section 3) as a function of the time step counter NC. Figure 
5-6 displays the net mass flow through the nozzle (i.e., the difference between 

the mass flows through the nozzle inlet and exit planes), normalized by the product 
of the stagnation density, stagnation sound speed, and the square of the nozzle 
throat half-height. The net mass flux should vanish at ateady state. The final 
mass flow is out of balance by over 1.2%, and it is apparent that the solution has 
not yet converged. The dimensionless thrust Fg shown in Fig. 5-7 also shows some 
variation throughout the course of the run. 


57 



GF ( 1 ) 

- 0.06 - 0.05 - 0.04 - 0.03 - 0.02 - 0.01 0.00 0.01 0.02 0.03 0.04 0.05 


NZLT2E FSMRCH- 1.000 RE=9.4+05 SMU-.050 LMRXBC= 1 L1BC= 1 



Figure 5-6 Convergence History of Net Mass Flow for 2-D Nozzle with Laminar 
Flow. Constant Time Stepsize At = 0.05. 
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Figure 5-8 Convergence History of Net Mass Flow for 2-D Nozzle with Laminar 
Flow. Variable Time Stepsize Increasing from Ax = 0.05 to 0.3. 
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Figure 5-9 Convergence History of Net Thrust for 2-D Nozzle with Laminar 
Flow. Variable Time Stepsize Increasing from At = 0.05 to 0.3. 
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Figures 5-8 and 5-9 display the corresponding results from a second 500 step run 
in which the stepsize Ax was increased automatically by the factor 1,1 after 
each step during which the maximum relative change | Aq/q | max was less than 1%. 

The stepsize increased gradually from its initial value Ax = 0.05 to a maximum value 
Ax = 0.3 during the first 250 steps, and remained at that value for the remaining 
250 steps. The largest stepsize Ax = 0.3 corresponds to a Courant number Co ( c ) = 
2600. The solution is essentially at steady state after 500 such increasing steps. 
The final net mass flow is out of balance (F^ = 0) by less than 0.001%, and the 
thrust is steady to five significant figures. Additional results from this 
converged case are given later in Section 6.1 and compared with the wind tunnel 
data. 


5.3 EFFECT OF ARTIFICIAL EXPLICIT SMOOTHING ON THE SOLUTION 

Earlier numerical experiments with the implicit algorith showed that the solution 
is not entirely insensitive to the explicit fourth-order smoothing terms which 
are incorporated to control short wavelength spatial oscillations (Ref. 1, Sections 

4.3 and 5.3). An array of numerical experiments has been performed to study 
parametrically the variation in the steady state solution when the magnitude of 
the smoothing coefficient is changed. 

A set of three smoothing terms is added explicitly to the R.H.S. of Eq. (2.16a). 
Each term represents a conservative fourth-order difference operator acting in 
one of the three coordinate directions £,n>C- According to a linear stability 
analysis (Ref. 1, Section 4.3), the coefficients of these terms are subject to 
the stability criterion 



(3-jAx 



< 1 


(5.8a) 


K 


i 



(5.8b) 


where x^, i = 1, 2, 3 represent the three directions £,ri)£ s (3 At denotes the 

i 

smoothing coefficient for each direction, and N = 2 or 3 for two or three- 
dimensional flow. In the numerical implementation, the product £ [l+ k - ] is 
taken to be the same for all directions. The stability bound on 1 the coefficient 
for each direction then is 
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3-jAx < 1/8N(1 +k 1 ) 


(5.9) 


We define a relative smoothing coefficient SMU as 

SMU = 3iAx/(3jAx) max (5.10) 

where (3jAx) max denotes the upper bound obtained by invoking the equality in 
(5.9). The quantity SMU is used as an input to the computer program to specify 
the smoothing coefficient (3-jAx) for each direction as a fraction of the maxi- 
mum stable value determined by inequality (5.9) 

3-jAx = SMU/8N(l+K-j ) , 0 < SMU < 1 (5.11) 

Note that, for a fixed value of SMU, the magnitude of the explicit smoothing 
terms relative to the actual spatially differentiated terms in the Navier-Stokes 
equations (2.1) depends upon the time stepsize Ax. The reason is that the 
smoothing terms are appended to the R.H.S. of Eq. (2.16a), in which the true 
spatially differentiated terms of the Navier-Stokes equations are multiplied by 
the factor Ax. This fact has been accounted for in Eq. (5.8) by defining the 
coefficient (3-jAx) of the smoothing term for each direction to include Ax as a 
factor. The magnitudes of the artificial smoothing terms relative to the actual 
terms in the Navier-Stokes equations then depends on the quantities 3-j, and would 
be independent of the time stepsize if each 3-j were taken as constant and inde- 
pendent of Ax. However, this would cause the stability criteria (5.8) or (5.9) 
to be violated for sufficiently large values of Ax. In terms of the quantity 
SMU, the stability criteria require that the coefficient 3-jAx satisfy Eq. (5.11). 
Since 3i itself governs the magnitude of the artificial smoothing terms relative 
to the actual Navier-Stokes terms, it follows from Eq. (5.11) that these relative 
magnitudes will be preserved under a change in Ax only when the ratio SMU/Ax is 
held fixed, rather than when SMU itself is held fixed. That is, the ratio SMU/Ax 
and not SMU itself governs the relative magnitudes of the smoothing terms and the 
true Navier-Stokes terms. 

We have performed an array of numerical experiments for the adiabatic flat plate 
laminar flow problem to determine the sensitivity of the steady-state solution 
to the magnitude of the smoothing strength SMU/Ax for a range of freestream Mach 
numbers 0.1 <_ <_ 3. The Reynolds number is 10 5 , the Prandtl number is unity, 

and viscosity is proportional to the temperature. 
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Table 5-2 displays the variation in the computed steady-state drag coefficient Cq 
with smoothing coefficient for various freestream Mach numbers. For each M^, the 
last column in the table gives the percentage by which the drag coefficient for a 
given value of SMU/At differs from the drag computed with SMU/At = 1. The results 
indicate that the solution is fairly insensitive to the magnitude of the smoothing 
coefficient, except at the lowest Mach number considered. For very low-speed flows, 
it appears that the accuracy of the solution is degraded significantly unless the 
ratio SMU/At is of the order of unity or less. 

The smoothing coefficient also was observed to have some effect on the rate of 
convergence to steady state at both the lowest and highest Mach numbers considered. 
Transient oscillations in C D were damped out somewhat more quickly with SMU/At = 10. 
In the supersonic case, = 3, such behavior is to be expected because the 
smoothing terms supply the dissipative mechanism that enables computation of the 
embedded shock wave induced by the displacement effect of the viscous boundary 
layer [see Ref. 1, Section 5 . 2 ]. 

5.4 EFFECT OF ARTIFICIAL IMPLICIT DISSIPATION ON STABILITY AND CONVERGENCE 

It has been noted in Section 5.2.2 above that, for the "parabolized" Navier-Stokes 
equations (2.1) wherein the viscous terms associated with the streamwise coordinate 
5 are neglected, the implicit algorithm tends to become unstable for large time 
steps At at which the streamwise Courant number Co(C) is much greater than the order 
of unity. We have performed some numerical experiments which show that the intro- 
duction of artificial implicit dissipation (see Section 2.4) extends the regime of 
stability to Cournat numbers Co(£) » 1 when the dissipation coefficient a is of 
the order of unity. This allows steady-state solutions to be attained in fewer 
time steps because the rate of convergence is generally greater the larger the 
time stepsize. 

The results of these experiments for adiabatic flat plate laminar flow at various 
freestream Mach numbers are summarized in Table 5-3. The table contains the 
results of two runs for each Mach number Mco. The first run employed no implicit 
dissipation (a = 0) and a stepsize At for which the streamwise Courant number Co(£) 
is of the order of unity. The second run used a dissipation coefficient a = 1 
and a variable step At, starting with the same initial flowfield and the same 
initial values of At and SMU (the explicit smoothing coefficient) as in the first 
run. The time step during the second run was increased automatically by the factor 
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Effect of Explicit Smoothing Coefficient on Steady State Drag Coefficient Adiabatic 
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Effect of Implicit Dissipation (a>0) on Stability and Convergence of Adiabatic Flat Plate* 
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1.1 following each step during which | Aq/q | max was ^ ess than unt ^ the 

time stepsize reached a pre-assigned maximum value Ax^, at which point the 
stepsize was held fixed for the remainder of the run. Whenever At was raised 
in this fashion, the explicit smoothing coefficient SMU also was raised by the 
same factor until it reached a maximum value SMU = 0.8. The latter bound was 
imposed to avoid violating the stability criterion for the smoothing terms 
j^Eq . (5.11) of Section 5.3]. The various columns in Table 5-3 list in order 
the freestream Mach number M*,; the initial and final stepsizes Ax^ and Ax-p; 
the streamwise Courant number Co.p(£) as estimated from Eq. (5.5) using Ax.p; 
the Courant number Co^t) based on Axf, on the maximum eigenvalue of 3h/3q, 
and on the minimum grid spacing in the wall normal direction; the final value 
of the explicit smoothing strength (SMU/At)^; the final drag coefficient Cpf 
at the last time step of the run; the time step index (NC)f at which the time 
stepsize reached its maximum value ATf, and the time step index (NC) S at 
which the drag coefficient first attained a value within +0.1% of Cpf and 
stayed in that range for the remainder of the run. This last entry in the 
table has been included merely to provide some quantitative measure of the 
convergence rate. Note that, for a given Mco, the final values of C D f differ 
slightly between the runs with fixed and variable stepsize Ax because the 
ratio (SMU/Ax) f could not be held fixed and still satisfy the stability 
criterion (5.11) in the variable-step run. 

For each freestream Mach number in the table, a run of fixed stepsize Ax for 
which Co(£) ~ 10 was attempted without implicit dissipation (a = 10) and was 
found to be unstable. In contrast, all of the variable-step runs with a = 1 
remained stable and converged rapidly for a maximum value of Co^(C) ~ 3000. 

The drag coefficient histories for these three runs are displayed in Figs. 5-10 
to 5-12. We haven't determined the maximum value of Co(£) that can be attained 
with a = 1 for various Mach numbers. However, an attempted variable-step run 
at 1^ = 0.1 with no pre-assigned bound on Ax eventually became unstable at a 
value Co(£) somewhat in excess of 4000. 

Implicit dissipation is potentially useful only for obtaining steady-state solutions. 
The artificial terms destroy the time-accuracy of the computation and lead to a 
physically unrealistic transient behavior. Neither does the implicit dissipation 
always ensure either stability or faster convergence to steady state. For 
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Figure 5-11 Flat Plate Drag Coefficient History Computed with Implicit 
Dissipation a = 1, (At) max = 100. FU = 0.8, Re*, = 10 5 . 


69 



GF [ 2 ) *10' 


o BL8B FSMRCH- 1.500 RE-1.0+05 SMU=.030 LMflXBC= 2 L1BC= 1 



Figure 5-12 Flat Plate Drag Coefficient History Computed with Implicit 
Dissipation a = 1, (Ar) max = 70. Moo = 1.5, Re^ = 10 . 
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example, we performed a run at Moo = 3 using a = 1 and a constant stepsize 
At = 0.01^Co(£) = 0.56J that failed to reach steady state in 200 steps; whereas, 
as indicated in line 5 of Table 5-3, the corresponding run with a = 0 converged 
adequately in the same number of steps. A similar run for Mm = 3 with an increasing 
stepsize became unstable at At = 0.02. This suggests that the implicit dissipation 
described in Section 2.4 may have a destabilizing effect for supersonic flows 
with embedded shock waves. In another example, an internal nozzle flow computation 
starting from very crude initial flowfield conditions proved unstable with a = 1 
for a small, fixed stepsize At | Co(£) * lj, but was stable with a = 0. In the 
light of the results shown in Table 5-3, these examples suggest that the implicit 
dissipation can enhance both stability and convergence for subsonic and transonic 
flows if very large time steps are employed, and if the initial conditions either 
are approximately correct, or are first regularized by a short preliminary run 
with a = 0 using shorter time steps. The applicability of implicit dissipation 
to supersonic flows with embedded shock waves requires further investigation. 

5.5 CONCLUSIONS 

The principal conclusions that have emerged from the numerical experiments described 
in the preceding subsections are summarized below. 

1. It has been demonstrated that the implicit numerical method is capable of 
computing viscous flows with boundary layer separation. 

2. Implicit boundary conditions are superior to explicit time-lagged boundary 
conditions in that they allow the use of larger time stepsize without 
numerical instability. For external flows, implicit freestream boundary 
conditions wherein the momentum components transverse to the boundary are 
computed implicitly from the transverse momentum equations are vastly 
superior to the simple method where all freestream conditions are imposed 
at the freestream boundary. The latter conditions retard convergence to 
steady state, and seriously degrade the accuracy of the computed flowfield. 

3. In general, larger time stepsizes give faster convergence to steady state. 
Although a linear stability analysis indicates the implicit algorithm to be 
unconditionally stable, experience indicates that the algorithm tends to 
become unstable at time stepsizes for which the streamwise Courant number 
exceeds the order of unity. Numerical experiments indicate that, within 
the stable range, the convergence rate is quadratic in the time stepsize. 
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4. The algorithm involves artificial explicit fourth order smoothing terms. 
Numerical experiments show that the accuracy of the solution is fairly 
insensitive to the magnitude of the smoothing coefficient, except for 
very low-speed subsonic flows. For such flows, the accuracy is degraded 
significantly unless the ratio of the smoothing coefficient to the time 
stepsize is of the order of unity or less. 

5. Numerical experiments indicate that, under favorable conditions, the 
introduction of artificial implicit dissipation extends the regime of 
numerical stability to streamwise Courant numbers far in excess of unity. 
This yield faster convergence to steady state when large time stepsizes are 
used (Co » 1). However, the implicit dissipation must be used with 
caution. It destroys the time-accuracy of the solution and cannot be 
employed in cases where the unsteady behavior of the flow is of interest. 
The present numerical experiments suggest that the implicit dissipation 

is most effective for subsonic and transonic flows. Further investigation 
is recommended to study its behavior for supersonic flow with embedded 
shock waves. 
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SECTION 6 


NOZZLE FLOWFIELD PREDICTIONS AND COMPARISON WITH EXPERIMENTAL DATA 


Laminar and turbulent flowfield computations have been performed for 
several three-dimensional nozzles that have been tested experimentally in the 
NASA Langley 16-foot Transonic Wind Tunnel. The first computation is for internal 
flow in a so-called "two-dimensional" converging-diverging nozzle with flat 
sidewalls and a rectangular cross-section. The second computation is for the 
combined internal and external flowfields of a circular, converging nozzle 
in a high subsonic external flow. The results of these computations are 
presented below. 

6.1 Internal Flow in a "Two-Dimensional" Converging-Diverging Nozzle 

6.1.1 Configuration, Operating Conditions, and Computational Grid 

The nozzle configuration and dimensions are shown in Fig. 6-1. The 
configuration is bilaterally symmetric with flat sidewalls and a rectangular 
cross section. The inlet section is of constant height. This is followed 
by a straight-walled converging section that is connected to a straight- 
walled diverging section by a circular arc that forms the geometric throat 
region. The exi t-to-throat area ratio is 1.0891 for a design exit Mach 
number of 1.35.' The corresponding design exit pressure and temperature from 
one-dimensional isentropic flow theory are p = 0.337 and T = 0.733, 
normalized by the stagnation chamber pressure and temperature, respectively. 

The operating conditions correspond to the design condition with a stagnation 
pressure of 1 atm. and a stagnation temperature of 295 Kelvin. The Reynolds 
number based on stagnation conditions and throat half-height is 9.3 x 10 5 as 
computed from the Sutherland viscosity law. The working fluid is air (y = 1.4, 
Pr = 0.72). 

For the flow computation, the origin of the Cartesian coordinate system 
is positioned at the geometric center of the throat. The x axis coincides 

with the intersection of the vertical and horizontal symmetry planes. 
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Figure 6-1 Configuration and Coordinates of Two-Dimensional Nozzle 



The latter coincide with the coordinate planes y = 0 and 2 = 0. Symmetry 
boundary conditions are applied at these planes, and the flow is computed in 
the quarter-space y>0, z>0. The inflow boundary is positioned at x = -4 in., 
and the outflow boundary is at the nozzle exit. The flow region between 
these planes, the symmetry planes, and the nozzle walls is covered by a 23 x 
10 x 15 grid in the x, y, and z directions, respectively. The streamwise (x) 
grid spacing is nonuniform with a relatively fine spacing near the throat, 
and the grid in each cross-sectional plane is exponentially stretched in both 
the y and z directions to resolve the wall boundary layers. The computational 
coordinate system thus has the form E, = 5(x), r\ = n(y), £ = c(z). The 
transverse grid spacings Ay, Az are several thousand times smaller at the 
walls than at the symmetry planes. Side and end views of the grid are 
displayed in Fig. 6-2a and 6-2b, and an orthographic projection of the 
three-dimensional grid is given in Fig, 6-2c. All dimensions shown are 
referred to the throat half-height. 

Laminar and turbulent flow computations have been performed for two- 
dimensional flow in the vertical plane of symmetry, and a laminar flow 
calculation has been performed on the three-dimensional grid shown in Fig. 6.2. 
All computations employ adiabatic wall boundary conditions. Initial conditions 
are obtained as follows. The inviscid core flow is obtained from one- 
dimensional isentropic flow theory for the nozzle area variation. These core 
flow conditions are applied over the central part of the grid in each cross- 
sectional plane assuming that the local velocity vector is oriented 

along the streamwise grid lines. Velocities on the remaining part of the 
grid are linearly interpolated in k and H to zero at the walls. Pressure is 
taken as uniform over each cross-sectional plane. Temperature in the nonuniform 
velocity region near the walls is obtained from velocity through the Crocco 
relation, and density follows from the equation of state. The resulting total 
pressure, total temperature, and transverse (v, w) distributions over the 
inflow plane are used as inflow boundary conditions for the flowfield compu- 
tation (see Section 2.3.1). The numerical results of the computations are 
presented below in dimensionless form. Dimensions are referred to the throat 
half-height, which is the reference length in terms of which the Reynolds 
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number is defined. Pressure, density, and temperature are referred to the 
corresponding stagnation conditions, and velocities are referred to the 
stagnation sound speed. 

6.1.2 Two- Dimensional Flow Results 

Laminar Flow 

The convergence history of this laminar flow computation has already 
been discussed in Section 5.2. Convergence was attained in 500 variable time 
steps. The first 100 of these were of constant size At = 0.05, the next 70 
were of increasing size, and the final 280 were constant at At = 0.3[Co(c) = 2600]. 
The computation took about 3 1/2 min. of CPU time on a CDC 7600 machine. The 
computed discharge coefficient is (^ = 0.9972, and the dimensionless thrust 
is 7.2612. The latter is referred to p 0 (c 0 L) 2 /2 where the subscript denotes 
stagnation conditions, c is the sound speed, and L is the reference length 
(the throat half-height). 

The nozzle wall pressure distribution is shown in Fig. 6-3, along with 
preliminary experimental data from the wind tunnel test (Ref. 10). The two 
sets of data represent pressure measurements on the upper and lower walls of 
the test nozzle. The computed pressure distribution has the same qualitative 
behavior as the data, including a slight recompression downstream of the 
geometric throat. The computation is in good quantitative agreement with 
the data in the subsonic and transonic region, but systematically under- 
predicts the data by about 5% in the supersonic region downstream of the 
throat. We shall see later that a turbulent flow computation agrees slightly 
better with the data in this region. 

The computed pressure distributions along the wall and the nozzle center- 
line displayed in Fig. 6-4 show a substantial difference in the subsonic and 
transonic region. 

Turbulent Flow 

The turbulent flow computation employed the converged laminar flowfield 
as initial conditions. The initial stepsize At = 0.05 increased by a factor 
of 6 over the first 75 steps and was held constant for 325 more steps. 
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Figure 6-3 Comparison of 2-D Laminar Flow Wall Pressure Distribution With 
Experimental Data of Ref. 10. Data Points: Circles, Upper 
Flap ; Triangles, Lower Flap. Computation: Solid Line. 
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although convergence was attained well before the end of the run. The 
computer time per step was about 13% greater than in the laminar flow 
calculation because of the extra computational labor involved in evaluating 
the turbulent viscosity. 

The computed steady-state discharge coefficient is = 0.9961, and 
the dimensionless thrust is 7.2922. The convergence history is displayed in 
Fig. 6-5 which shows the dimensionless net mass flux and thrust as a function 
of the time step index NC. 

The computed wall pressure distribution displayed in Fig. 6-6 is in 
somewhat better agreement with the experimental data than the laminar flow 
results shown in Fig. 6-3. The computed wall and centerline pressure distri- 
bution are shown in Fig. 6-7. The vertical profiles of velocity and pressure 
across the nozzle throat are given in Fig. 6-8. Similar profiles across the 
nozzle exit plane (the outflow boundary) are given in Fig. 6-9. Figures 
6-7 to 6-9 show that substantial flowfield gradients exist across the nozzle 
in the neighborhood of the throat. In the supersonic region downstream of 
the throat, however, the pressure becomes uniform across the nozzle* The 
velocity is also spatially uniform, except in the near-wall region occupied 
by the turbulent boundary layer. 
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Figure 6-5a Convergence History for 2-D Turbulent Flow Computation Showing 
Net Mass Flow Rate Versus Time Step Index NC. 
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Figure 6-6 Comparison of 2-D Turbulent Flow Wall Pressure Distribution With 
Experimental Data of Ref. 10. Data Points: Circles, Upper Flapj 

Triangles, Lower Flap. Computation: Solid Line. 
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Figure 6-7 Computed Wall and Centerline Pressure Distributions for 2-D 

Turbulent Flow. Solid Line: Wall. Dashed Line: Centerline. 
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Figure 6-8a Velocity Profiles Across Throat for 2-D Turbulent Flow 
Computation 
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Figure 6-8b Pressure Profile Across Throat for 2-D Turbulent Flow Computation 
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Figure 6-9a Velocity Profiles Across Exit Plane for 2-D Turbulent Flow 
Computation 
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Figure 6-9b Pressure Profile Across Exit Plane for 2-D Turbulent Flow 
Computation 
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6.1.3 Three-Dimensional Flow Results 

A three-dimensional laminar flow computation was performed on the 
grid shown in Fig. 6-2. The computation was run for 500 time steps, by which 
time the computed pressure field was essentially converged and was in agreement 
with the two-dimensional flow solution at the vertical plane of symmetry. The 
maximum Courant number (based on the minimum transverse mesh spacing) varied 
from 94 to 516 over the course of the run, which took approximately one hour 
of CDC 7600 computer time. The computed net thrust and discharge coefficient 
are 7.2197 and 0.9774, respectively. Both are somewhat lower than the two- 
dimensional flow predictions, which do not account for the sidewall. 

Figure 6-10 shows a comparison of the computed wall pressures with 
experimental data. The pressure distribution along the upper nozzle wall at 
the vertical plane of symmetry is given in Fig. 6-10a. The corresponding 
pressure distribution along the sidewall at the horizontal plane of symmetry 
is displayed in Fig. 6-10b, As for the two-dimensional results, the compu- 
tation tends to underpredict the data downstream of the throat, but agrees 
closely with the data in the subsonic and transonic regions upstream of the 
throat. 

The computed flowfield is smooth, regular, and shows no evidence of 
boundary layer separation over most of the nozzle interior. However, the 
flowfield contains a strong secondary flow in the neighborhood of the axial 
corner where the upper wall and sidewall intersect. This secondary flow is 
localized in the throat region -0.88 £ x < 0.6 (see Fig. 6-2a), and displays 
a flow reversal. Near the upstream and downstream ends of the described x 
interval, the velocity profiles near the axial corner display a slight flow 
reversal that is reminiscent of the reverse flow profile in a separated 
boundary layer. However, there is no flow reversal over the bulk of the 
upper wall or sidewall except in the immediate neighborhood of the corner. 

The peak reverse velocity exists at an axial station just upstream of the 
throat (x = -0.186), and is about 100% greater in magnitude than the general 
streamwise velocity V c =s0.8 along the nozzle centerline at the same x station. 
The secondary flow thus is more properly characterized as a reverse jet rather 
than as a region of separated flow. The jet-like structure is readily apparent 
in Figs. 6-lla, b, which show the computed velocity profiles in the cross- 
sectional plane x = -0.186. Figs. 6-lla, b display the velocity profile along a 
vertical grid line Y = const, that passes approximately through the peak 
velocity point of the jet. Figs .. 6rll c, d display a similar profile along a 
horizontal grid line Z = const, that passes approximately through the peak. 
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Figure 6-10a Computed Pressure Distribution Along Upper Nozzle Wall 
at Vertical Symmetry Plane Compared with Experimental 
Data 
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Figure 6-10b Computed Pressure Distribution Along Nozzle Sidewall at 
Horizontal Symmetry Plane Compared with Experimental Data 
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Figure 6-llb Magnified View of Near-Wall Portion of Velocity 
Profile Shown in Figure 6-lla 











Figure 6-1 Id Magnified View of Near-Wall Portion of Velocity Profile 
Shown in Figure 6-llc 




6.2 INTERNAL AND EXTERNAL FLOWFIELD OF A CIRCULAR NOZZLE 

6.2.1 Configuration, Operating Conditions, and Computational Grid 

The configuration of this ax i symmetric nozzle is displayed in Fig. 6-12. 

The interior wall has an inlet section of constant radius, a convergent section, 
and an exit section of constant radius. The external wall consists of an initial 
section of constant radius, followed by a circular arc boattail. The internal 
and external surfaces do not meet at a sharp trailing edge at the exit. Rather, 
the trailing edge is squared off by a vertical cut. Because the numerical 
algorithm employs a Cartesian base coordinate system, axi symmetric flows must 
be treated as any other three-dimensional flow; that is by using a full three- 
dimensional grid. In the present case, the flow region about the nozzle is 
bounded by inflow and outflow planes normal to the symmetry axis, and by a 
cylindrical outer boundary (the freestream boundary) that is concentric with the 
nozzle and is located 5 cm away from the cylindrical initial section of the 
outer surface of the nozzle. The outflow boundary is placed about one exit 
diameter downstream of the nozzle exit so as to include the near-field exhaust 
jet and the wake of the nozzle wall in the computation. The first quadrant y>0, 
z>0 of the region between the inflow and outflow planes and within the outer 
cylindrical boundary is mapped onto a rectangular computational parallelepiped 
as shown in Fig. 6-13. The flow computation is carried out on a uniform 
28x5x28 grid in this computational space. Note that the mapping is singular, 
inasmuch as the axis of symmetry maps onto the surface c=0 in the computational 
space. 

The image of this uniform grid in physical space consists of five meridional 
planes (including the symmetry planes y=0 and z=0) equally-spaced at intervals of 
22.5 degrees; 28 cross-sectional planes along the x direction, of which the last 
five are downstream of the nozzle exit; 15 grid points distributed exponentially 
across the interior of the nozzle in the radial direction, and another 13 points 
distributed exponentially between the outer surface of the nozzle wall and the 
cylindrical outer boundary. Side and end views of the physical grid are displayed 
in Figs. 6-14a and 6-14b. 

The wind tunnel test conditions are such that both the outer and internal 
flows are turbulent. The conditions are as follows (Ref. 11): Air is the medium 
for both the internal and external flows (y = 1.4, Pr = 0.72). The internal flow 
stagnation pressure and temperature are 1.32 atm. and 300 Kelvin, respectively. 
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Figure 6-12 Configuration and Coordinates for Circular Nozzle 
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Figure 6-13 Illustration of Computational Space for Circular Nozzle 
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Figure 6-14b End View of Computational Grid for Circular Nozzle 


102 


The external flow has a freestream Mach number Mo = 0.8, a stagnation temperature 
of 275 Kelvin, and the freestream static pressure is half the internal flow 
stagnation pressure. The Reynolds number based on internal flow stagnation 
chamber conditions and on the nozzle exit radius is 1.1 x 10^, using the 
Sutherland viscosity law. 

The boundary conditions for the computation are as follows: Adiabatic 

wall boundary conditions are employed. The transverse velocity components 
are assumed to be small over the entire inflow boundary (v = w = 0), and the 
total temperature over the external flow portion of the inflow boundary is taken 
as uniform and equal to the freestream stagnation temperature. Static pressure 
surveys taken during the test showed no significant departure from the freestream 
pressure over the external flow part of the inflow plane (Ref. 12). This 
information, together with the measured external boundary layer velocity distri- 
bution has been used to calculate the remaining flow variables, which are held 
fixed over the external flow portion of the inflow plane during the flowfield 
computation. 

On the portion of the inflow plane interior to the nozzle, probe surveys 
lacked the resolution to give information any more detailed than an estimate 
on the boundary layer thickness on the interior wall surface (Ref. 12). 
Consequently, the total temperature has been assumed uniform and equal to the 
measured stagnation temperature. A one-seventh power boundary layer velocity 
profile has been used together with the static pressure and core velocity 
predicted by one-dimensional isentropic flow theory to estimate the total pressure 
distribution across the part of the inflow plane interior to the nozzle. The 
described total pressure and total temperature are used as boundary conditions 
for the implicit flowfield computation (see Section 2.3.1). 

The initial conditions for the region between the inflow plane and the 
nozzle exit plane are obtained as follows. Flow conditions in the interior of 
the nozzle are estimated in the same fashion as described in Section 5.1 for 
the cosine-shaped nozzle. Outside the nozzle, the transverse velocity components 
are set to zero (v = w = 0), and both the static pressure and the total tempera- 
ture are assumed uniform and equal to their freestream values. Along each 
streamwise mesh line n = const., c = const., the total pressure is assumed 
constant and equal to its value at the inflow plane. The velocity u and static 
temperature T then are computed from the known total pressure, total temperature, 
and static pressure, and density follows from the equation of state. The 
initial conditions at grid points in the cross sectional planes downstream of 
the nozzle exit plane simply are set equal to the conditions at the corresponding 
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points in the exit plane itself. 

Numerical results for turbulent flow are presented below in dimensionless 
form. The non-dimensionalization is as described at the end of subsection 6.1.1. 

6.2.2 Numerical Results 

The three-dimensional turbulent flow computation attained convergence 
in 500 time steps at maximum Courant numbers (based on the minimum transverse 
grid spacing) that ranged from 1400 to 6800. The computation took slightly less 
than one hour of CDC 7600 computer time, and yielded a final steady state value 
of 3.428 for the net thrust. 

Figure 6-15 displays the computed static pressure distribution along the 
interior surface of the nozzle wall and along the continuation of the wall grid 
line into the downstream mixing region. A similar distribution of pressure along 
the exterior surface of the nozzle wall is given in Fig. 6-16. Pressure co- 
efficients derived from the latter data are shown in Fig. 6-17, along with wind 
tunnel experimental data from pressure taps located along the boattail (Ref. 11). 
The qualitative character of the theoretical curve matches that of the data, but 
the two are in close quantitative agreement only near the aft end of the boattail. 
The theoretical prediction of the crossover point from negative to positive values 
of Cp lies somewhat downstream of the experimental crossover point, and the 
minimum value of Cp that occurs near the forward end of the boattail is lower in 
the experimental data. A similar quantitative discrepancy with data on the 
forward part of the boattail was found in an axisymmetric flow solution by inves- 
tigators at the Langley Research Center (according to Lawrence Putnam). In that 
solution, the jet mixing region was modeled as a cylindrical sting, and only the 
external flow was computed. Numerical experimentation showed that the theoreti- 
cal prediction could be brought into agreement with the data by moving the upper 
computational boundary (freestream boundary) farther away from the nozzle outer 
surface. In the present computation, that boundary is only about 1.3 nozzle 
exit radii above the forward end of the boattail, whereas the corresponding dis- 
tance is 2.3 exit radii at the aft end of the boattail where the theory and 
experiment agree closely. This suggests that the disagreement in the forward 
region could be eliminated by shifting the freestream boundary outward a dis- 
tance of about one exit radius. 
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The computed flowfield shows that the external boundary layer remains 
attached over the most of the boattail, although there is a small region of 
separated flow near the aft end; specifically, there is a slight flow reversal 
in the velocity profiles at the last two streamwise grid points on the boattail, 
as shown in Fig. 6-18. It should be noted that an earlier laminar flow compu- 
tation for the same geometry and operating conditions displayed massive 
separation over most of the boattail. 

In addition to the boattail pressure taps, the wind tunnel test data 
included total pressure probe surveys across several axial stations in the 
downstream jet mixing region (Ref. 13). Figure 6-19 shows a comparison of these 
data with the computed total pressure distributions along the radial grid lines 
that lie nearest the probe survey stations, normalized by the computed values 
at the outermost interior grid point. One can see that at each station the 
computed shear layer in the near wake of the nozzle wall is somewhat thinner 
than that measured experimentally. This suggests that the wake turbulence 
model described in Section 4 may tend to under-predict the magnitude of the 
turbulent eddy viscosity in such regions. 

One can also see from Fig. 6-19 that the computation apparently over- 
predicts the total pressure in the inner core of the jet. This is a result of 
the coarse grid employed in the region near the symmetry axis. It is anticipated 
that the agreement between predicted and experimental results would be improved 
substantially with the same number of grid points by using a less highly 
stretched radial grid (the grid spacing at the symmetry axis is a factor of 
5000 greater at the axis than at the interior nozzle wall in the present 
computation), and by shifting the freestream computational boundary outward by 
a distance of one nozzle exit radius. 
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X(J) 

Figure 6-15 Computed Pressure Distribution Along Interior of Nozzle 
Wall and Along Continuation of Wall Grid Line into 
Downstream Jet Mixing Region 
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Figure 6-16 Pressure Distribution Along Exterior of Nozzle Wall and 
Along Continuation of Wall Grid Line into Downstream 
Jet Mixing Region 
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Figure 6-17 Pressure Coefficients Derived From Pressure Distribution of Figure 6-16 
Compared with Wind Tunnel Data 
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Figure 6-18 Velocity Profiles Along a Vertical Grid Line at 
Nozzle Exit Plane, x = 2.0 
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Figure 6-19b Comparison Between Experimental Total Pressure 
Probe Survey at x=3.128 with Normalized Total 
Pressure Profile Along the Nearest Vertical 
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